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This  work  explores  several  highly  energetic  molecular  systems  that  might  be 

useful  as  novel  fuels  or  propellants.  These  systems  include  04,  N603,  and  a  series  of 

azacubanes  (C8-„NnH8_n)  along  with  their  their  nitro-substituted  derivatives 

(C8-nN„(N02)8-n)-  For  O4,  a  cooperative  experimental-theoretical  study  was 

performed  and  a  novel,  long-lasting,  van  der  Waals  complex  was  spectroscopically 

observed  and  identified.  However,  another  isomer  of  O4,  a  covalently-bound  D3/l 

form,  seems  promising  as  a  highly  energetic  metastable  species  and  merits  further 

study.  The  N603  system  was  studied  in  terms  of  a  series  of  nitrogen  rings  stablized 

by  coordinate-covalent  bonds  with  oxygen.  Out  of  all  of  the  nitrogen-ring  molecules 

considered,  N603  had  the  best  separation  of  adjacent  atomic  charges.  Consequently, 

it  also  had  the  most  kinetic  stability  with  a  62.4  kcal  mol-1  activation  energy 

toward  unimolecular  dissociation.  The  series  of  azacubanes  were  investigated  as 

derivatives  of  the  cubane  molecule,  an  extremely  dense,  shock-insensitive  explosive. 

The  computed  heats  of  formation  for  all  of  the  azacubanes  were  larger  than 

cubane's,  however,  none  of  the  azacubanes  have  been  synthesized. 


To  facilitate  predictive  calculations  on  other  highly  energetic  and  larger 
molecules,  several  improvements  for  the  virtual  orbital  space  were  explored.  These 
include  virtual  orbitals  generated  from  a  potential  of  a  reduced  number  of  electrons, 
virtuals  generated  from  density  functional  potentials,  as  well  as  the  use  of 
approximate  natural  orbitals  for  the  virtual  space.  Numerical  results  show  that 
natural  orbitals  are  one  of  the  better  sets  and  recover  most  of  the  correlation  energy 
in  a  much  smaller  space. 
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CHAPTER  1 
INTRODUCTION 

Electronic  structure  theory  has  emerged  as  a  powerful  tool  to  elucidate 
molecular  structures,  energetics,  and  properties.  Developments  in  the  field,  coupled 
with  improvements  in  computer  technologies,  have  made  it  fairly  routine  to 
characterize  molecular  systems  with  five  or  fewer  atoms.  Typical  quantities 
predicted  are  molecular  potential  energy  surfaces  including  minima  and  extrema, 
energy  differences  with  their  associated  transition  states  and  activation  barriers,  and 
normal  harmonic  vibrational  frequencies  with  their  infrared  and  Raman  intensities. 

The  coupled-cluster  (CC)  framework,  with  an  exponential  ansatz,  is  often  the 
method  of  choice  for  solutions  to  the  electronic  Schrodinger  equation,  as  well  as  the 
prediction  of  properties  [1-6].  With  the  demise  of  configuration-interaction  (CI) 
methods  in  the  1990s,  years  of  chemical  literature  show  that  CC  methods  are  both 
accurate  and  routinely  applicable.  Current  research  efforts  that  add  quadruple  [7], 
pentuple  [8],  and  higher  excitations  [9]  (as  well  as  the  efficient  construction  of 
spin-adapted  wave  functions  [10-12])  will  allow  CC  theory  to  describe  even  more 
chemical  processes. 

Coupled-cluster  methods  have  made  great  progress  in  the  field  of  high 
energy-density  materials  (HEDMs)  [13,14].  Potential  HEDMs  are  characterized  by 
highly  energetic  systems  that  are  metastable.  It  is  also  desirable,  but  not  required, 
that  HEDMs  be  environmentally  friendly.  One  example  is  liquid  hydrogen  and 
oxygen  which  is  used  as  a  propellant  in  the  space  shuttle  [15].  Other  HEDMs 
include  monopropellants  that  are  not  combined  with  oxygen  and  that  produce 
thrust  via  unimolecular  dissociation  to  gaseous  products. 


Coupled-cluster  methods  have  many  strengths  that  make  them  extremely 
useful  for  the  study  of  HEDMs.  The  first  is  their  ability  to  survey  many  different 
types  of  systems.  Provided  that  a  one-particle  Gaussian  basis  set  has  been 
developed  for  the  elements  of  interest  and  that  relativistic  effects  are  small  or 
included  by  pseudo  potentials,  CC  methods  should  be  applicable.  In  fact,  a 
theoretical  approach  facilitates  the  construction  of  HEDMs  by  design.  For  instance, 
atoms  can  be  easily  changed  within  the  input  and  substituent  effects  deduced. 
Furthermore,  theoretical  studies  of  environmentally  unfriendly  fuels,  propellants  and 
explosives  (i.e.,  those  that  produce  metals  or  toxic  products)  are  safer  and  can  be 
less  costly  than  physical  experiments.  Also,  CC  methods  can  be  highly  accurate  in 
terms  of  predicted  bond  lengths,  angles  and  normal  harmonic  frequencies,  and  thus 
predictive.  Often  potential  HEDMs  exploit  novel  bonding  patterns,  or  rather 
bonding  patterns  that  are  not  yet  known,  so  highly  accurate  methods  are  needed. 
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Figure  1-1.  The  computed  CCSD(T)/aug-cc-pVTZ  structure  of  tetrahedral  N4 
(tetrazete). 


One  example  of  a  molecule  envisioned,  characterized  [16-30],  and  probably 
spectroscopically  observed  [31]  within  the  HEDM  effort,  is  the  tetrahedral  N4 
molecule,  shown  in  Figure  1-1.  Theoretical  calculations  using  CC  methods  show 
that  the  dissociation  energy  for  N4  into  2  N2  molecules  is  183  kcal  mol-1  [23]  and 
the  activation  energy  for  this  process  is  61  kcal  mol-1  [27].  Furthermore,  if 
spin-forbidden  radiationless  decay  is  considered,  the  barrier  is  still  28.2  kcal 
mol-1  [25].  Armed  with  this  information,  several  experimental  groups  [31,32] 
attempted  to  prepare  and  spectroscopically  observe  N4.  Because  of  its  high 
symmetry  (Td),  one  of  difficulties  is  that  only  one  vibrational  mode  of  t2  symmetry 


is  weakly  active  in  the  infrared.  In  fact,  a  weak  transition  was  observed  at  936.7 
cm-1  which  compares  very  well  with  the  CCSD(T)/aug-cc-pVTZ  computed  value  of 
936  cm-1  [31].  However,  the  observed  and  computed  frequencies  for  the  isotopically 
substituted  15N4  were  somewhat  less  consistent  with  values  of  900.0  cm-1  and  904.4 
cm-1,  respectively.  Although  three  vibrations  (e,  t2,  and  ai)  are  active  in  its  Raman 
spectra,  a  higher  concentration  (  >  80  ppm)  in  solid  N2  is  required  for  detection  [32]. 
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Figure  1-2.  The  computed  CCSD(T)/aug-cc-pVTZ  structure  of  pentagonal  N5 
(pentazole). 

Another  molecule  of  interest  in  the  HEDM  community  is  the  pentazole  ring 
shown  in  Figure  1-2.  Although  arylpentazoles  were  prepared  over  forty  years 
ago  [33-36],  the  Ng"  moiety  has  not  yet  been  isolated.  While  it  is  certainly  true  that 
strong  electron  withdrawing  groups  would  stabilize  the  five-membered  nitrogen  ring, 
a  pure  nitrogen  anion  with  a  large  ionization  potential  is  vitally  important  for  the 
synthesis  of  larger  nitrogen  molecules.  Just  recently,  Christe  et  al.  reported  the 
synthesis  and  characterization  of  the  highly  energetic  N^  molecule  [37-39]. 
Subsequent  work  explored  the  prospects  of  combining  N^"  with  the  azide  anion 
(N3  ).  These  attempts  have  not  been  successful  and  a  theoretical  study  of  the 
possible  N8  products  showed  little  promise  [40,41]. 

With  a  much  larger  ionization  potential  of  5.3  eV  compared  to  that  for  azide  of 
2.5  eV,  the  pentazole  anion  offers  a  more  promising  route  toward  novel  nitrogen 
molecules.  First,  the  pentazole  anion  is  completely  aromatic  while  the  related 
neutral  structure  is  either  weakly  bound  or  completely  dissociative  [42].  Second, 
pentazole's  activation  barrier  for  unimolecular  dissociation  is  reasonably  high  (27.4 
kcal  mol"1  at  the  CCSD(T)  level  of  theory  [43-47])  and  slightly  larger  than  the  19.9 


kcal  mol-1  activation  barrier  for  hydrogen  pentazole  [46,48].  This  suggests  that 
pentazole  would  be  kinetically  stable  if  it  could  be  prepared.  Third,  pentazole  might 
exist  bound  to  a  metal  atom  and  a  recent  study  showed  that  in  plane  bidentate 
complexes  are  more  strongly  bound  than  metallocene-like  or  sandwich 
complexes  [49].  Interesting,  Gagliardi  and  Pyykko  [50]  have  proposed  a  larger  N?- 
ring  with  10  it  electons.  Such  a  system,  would  be  able  to  donate  electrons  from  the 
5  shell  to  the  metal. 

If  CC  methods  are  to  provide  accurate,  predictive  and  computationally 
affordable  structures,  energetics,  and  properties,  a  significant  limitation  is  the 
dimension  of  the  molecular  orbital  (MO)  basis. a    For  example,  the  number  of 
floating-point  operations  per  closed-shell  iteration  of  the  coupled-cluster  method 
with  inclusion  of  connected  single  and  double  cluster  operators  (CCSD)  is 

^  {\n^Nvir  +  KJ&)  (1-1) 

where  h  is  the  number  of  symmetry  operations,  nocc  is  the  number  of  occupied 
orbitals,  and  Nvir  is  the  number  of  virtual  orbitals.  Similarly,  the  coupled-cluster 
method  with  inclusion  of  connected  single  and  double  cluster  operators  augmented 
by  a  noniterative  inclusion  of  triple  excitations,  CCSD(T),  scales  as 

4^  KMr  +  <cMr)  (1-2) 

Given  such  a  power  dependence  (which  is  sometimes  referred  to  as  the  exponential 
scaling  wall  [53]),  it  is  easy  to  see  that  applications  quickly  become  prohibitively 
expensive  with  increasing  molecular  size. 


a  Often,  the  dimension  of  the  MO  basis  is  the  same  as  the  primitive  atomic  orbital 
(AO)  basis.  The  need  for  large  AO  basis  sets  with  high  angular  momenta  is  one  of 
the  consequences  of  the  electronic  Coulomb  cusp  condition  [51,52]. 


Conversely,  the  power  dependence  as  a  function  of  the  number  of  orbitals  also 
displays  the  potential  savings  if  a  subset  of  the  full  space  is  used.  In  this 
dissertation,  I  show  that  frozen  natural  orbitals  (FNOs)  of  a  reduced  dimension  can 
be  generated  such  that 

N'viT  =  xNvir        where        x  <  1  (1-3) 

Generation  of  FNOs  in  the  full  space  is  approximate  (i.e.,  done  at  a  lower  level  of 
theory),  and  hence  does  not  involve  a  rate  limiting  step.  Although  FNOs  are  of  a 
smaller  dimension  than  typical  canonical  Hartree-Fock  virtual  orbitals,  they  still 
recover  nearly  all  of  the  correlation  energy.  This  translates  to  a  savings  of  x4  per 
iteration  in  the  number  of  operations  for  CCSD  and  CCSD(T),  and  even  larger 
savings  for  higher-level  calculations.  Comparisons  among  SCF  virtual  orbitals, 
those  obtained  from  effective  Fn_i  and  Vn_a  potentials  where  a  is  the  number  of 
valence  electrons,  those  using  density-functional  theory  (DFT)  potentials,  and 
FNOs  will  be  presented. 


CHAPTER  2 
UNRAVELING  THE  MYSTERIES  OF  METASTABLE  OJa 

Interest  in  tetraoxygen  molecules  dates  to  a  1924  paper  by  G.  N.  Lewis  [54]; 

since  then,  studies  of  O2  dimers  have  enjoyed  a  long  history  of  investigation  [55-57]. 

Theoretical  studies  of  covalently  bound  O4  species  began  with  Adamantides'  1980 

prediction  [58]  of  a  bound  cyclic  (D2d)  form.  This  stimulated  considerable  further 

theoretical  effort  as  this  cyclic  O4,  nearly  4  eV  higher  in  energy  than  two  O2 

molecules,  appeared  to  be  a  promising  candidate  for  a  high  energy  density 

material  [59-61].  Subsequent  theoretical  studies  also  identified  a  Dsh  form 

analogous  to  S03  at  a  somewhat  higher  energy  [62,63].  Although  experimentalists 

have  long  observed  evidence  of  van  der  Waals'  complexes  of  ground  state  O2 

molecules,  no  evidence  has  been  found  supporting  the  theoretical  predictions  of 

covalent  04  species.  In  a  recent  report  from  the  Suits'  laboratory  [64],  1+1  resonant 

photoionization  spectra  were  reported  for  an  energetic,  metastable  O4  species 

produced  in  a  DC  discharge  [64].  Intense  spectra  were  observed  throughout  the 

region  from  280  to  325  nm,  implying  an  initial  form  of  tetraoxygen  containing  at 

least  4  eV  internal  energy  relative  to  02  +  O2,  as  well  as  the  existence  of  a  higher 

excited  state  through  which  the  ionization  takes  place.  In  the  absence  of  plausible 

alternatives  accounting  for  all  the  observations,  an  energetic  covalent  O4  species  was 

considered  the  most  likely  candidate.  They  presented  rotationally  resolved 


a  This  chapter  consists  of  parts  of  an  article  reprinted  with  permission  from 
Darcy  S.  Peterka,  Musahid  Ahmed,  Arthur  Suits,  Kenneth  J.  Wilson,  Anatoli 
Korkin,  Marcel  Nooijen,  and  Rodney  J.  Bartlett,  Journal  of  Chemical  Physics, 
volume  110,  pages  6095-6098.  ©1999  American  Institute  of  Physics. 
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Figure  2-1.  Various  possible  isomers  of  O4. 

photoionization  spectra  and  photoelectron  spectra.  Here,  we  present  the  theory 
component  that  provides  compelling  indirect  evidence  pointing  to  the  identity  of 
this  species  as  a  novel  complex,  involving  one  ground  state  02  molecule  and  one  in 
the  metastable  c  (*£")  state,  but  not  a  covalently  bonded  04. 

2.1     Experimental  Discussion 
The  experiment  [64, 65]  consists  of  passing  a  pulsed  molecular  beam  of  oxygen 
through  electrodes  held  at  ground  and  ±3  to  5  kV  so  that  the  discharge  occurs  in 
the  collision  region  of  the  beam.  The  discharge  is  positively  biased  when  detecting 
ions,  and  negatively  biased  when  detecting  electrons  to  inhibit  interference  from 
corresponding  species  in  the  beam.  The  molecular  beam  is  skimmed  before  entering 
a  main  chamber  wherein  it  is  crossed  by  an  unfocused  (for  wavelength  scans)  or 
loosely  focused  (for  photoelectron  spectra)  Nd-YAG  pumped  dye  laser  doubled  to 
yield  light  tunable  around  300  nm  with  a  linewidth  on  the  order  of  0.08  cm-1.  The 
laser  and  molecular  beams  cross  on  the  axis  of  a  time-of-flight  mass  spectrometer 
with  velocity  map  imaging  [66]  (VELMI)  detector,  allowing  for  several  different 
kinds  of  experiments  to  be  performed.  Mass-selected  resonant  ionization  scans  are 
effected  by  recording  the  total  mass-selected  ion  yield  as  a  function  of  laser 
wavelength.  Total  photoelectron  signals  are  similarly  obtained  by  reversing  the 
potentials  and  recording  the  integrated  electron  signal  striking  the  detector  as  a 
function  of  laser  wavelength.  Photoelectron  images  are  recorded  on  the  resonant 
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Figure  2-2.  Raw  0|  photoion  yield  and  total  photoelectron  yield  spectra. 

Expanded  region  of  the  electron  spectrum  is  shown  in  the  inset. 

lines  using  the  VELMI  technique,  calibrated  with  Ar*  ionization,  and  converted  to 
electron  kinetic  energy  using  established  techniques. 

Photoionization  spectra  recorded  for  m/e=64,  Oj  from  Suits  et  al.  [67]  are 
shown  for  the  region  from  302  to  325  nm  in  Figure  2-2.  Total  photoelectron  yield 
signals  are  also  shown;  these  signals  are  recorded  under  conditions  in  which  the 
discharge  bias  is  reversed  and  the  total  discharge  power  in  the  electron  case  is  lower 
(2.5  Watts  vs.  6  Watts).  The  experimentalists  ascribe  the  differences  in  the  ion  and 
electron  spectra  to  inherent  noise  in  the  higher-power  ion  scans  and  the 
correspondingly  higher  temperature  of  the  ion  scans,  giving  more  'hot  band' 
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Figure  2-3.  Relative  energies  (CCSD/TZ2P)  with  respect  to  02  X  (3Ej)  +  02  X 
(3E~)  along  different  ionization  paths. 


contributions.  Under  the  conditions  of  the  experiment,  virtually  no  ions  are 
observed  other  than  0|.  Inset  in  Figure  2-2  is  an  expanded  view  of  the  electron 
yield  spectrum  in  the  long  wavelength  region  near  323  nm.  Clearly  resolved 
rotational  spectra  are  observed,  with  line  spacings  on  the  order  of  3.2-3.6  cm"1. 
Similar  rotational  structure  is  also  apparent  on  some  lines  in  the  306  nm  region, 
although  not  as  pronounced. 
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2.2     Theoretical  Discussion 

If  one  of  the  covalent  species  is  responsible  for  these  spectra,  then  two  critical 
issues  are  1)  accounting  for  the  observed  ionization  potential  of  ~8  eV  from  the 
metastable  state  or  likely  12  eV  or  so  from  two  ground  state  O2  molecules,  and  2) 
finding  a  bound  excited  state  ~4  eV  above  the  metastable  species.  To  this  end,  we 
have  performed  accurate  coupled-cluster  (CC)  calculations  with  the  ACES  II 
program  system  [68].  We  use  a  TZ2P  basis  [69]  of  Cartesian  Gaussians  contracted 
as  (Ils6p3d)/[5s3p2d]  except  as  indicated.  Table  2-1  presents  computed  CCSD  and 
CCSD(T)  [1,2]  energies  relative  to  two  ground  state  O2  molecules  for  several  states 
of  interest.  At  the  CCSD(T)  level,  we  find  two  covalent  forms:  the  cyclic  (D2</)  at 
4.16  eV  and  the  pinwheel  (D3/l)  at  5.07  eV.  We  also  obtained  detailed  structures 
and  vibrational  frequencies  (shown  in  Table  2-2)  for  these  two  species  and 
IP-EOM-CCSD  [1,2]  vertical  ionization  potentials  for  the  two  covalent  forms.  The 
lowest  IPs  occur  at  10.98  for  the  cyclic  and  12.47  for  the  pinwheel  structures. 
However,  this  is  much  higher  than  the  energy  of  two  photons  at  300  nm  (~8  eV), 
thus  outside  the  range  of  the  experiment.  Also  shown  in  Table  2-1  are  the  energies 
of  the  ionized  van  der  Waals  complexes:  our  result  of  11.57  eV  is  in  excellent 
agreement  with  the  11.67  eV  CASSCF  value  of  Lindh  and  Barnes  [70].  Also  shown 
in  Table  2-1  are  rotational  constants  for  these  species.  These  are  not  consistent 
with  the  rotationally  resolved  spectra  in  the  inset  in  Figure  2-2.  Even  accounting 
for  the  nuclear  spin  symmetry  for  the  Dzh  species,  a  line  spacing  on  the  order  of  6Be 
or  1.2  cm-1  is  the  largest  expected. 

The  other  question  pertains  to  the  existence  of  a  bound  excited  state  of  the 
covalent  species  that  is  4  eV  above  the  metastable  species  and  whose  vibrational 
frequencies  were  previously  estimated  [64].  A  STEOM-CC  [71]  calculation  in  the 
POL1  basis  at  the  geometry  of  the  D2d  ground  state  shows  a  weakly  allowed  E  state 
at  7.10  eV  and  a  dipole  forbidden  A2  state  at  8.54  eV  and  three  others  weakly 
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Figure  2-4.  Relevant  potential  curves  adapted  from  Reference  [73],  from  calculations 
of  Reference  [74].  The  Rydberg  and  ion  curves,  duplicated  and  offset 
-0.45  eV  (the  energy  of  the  02-Oj  bond)  are  shown  as  dashed  lines. 

allowed,  between  9.2  and  10.0  eV.  For  the  D3h  form,  the  first  state  is  a  forbidden  Ax 
that  occurs  at  7.24  eV  with  a  strong  E'  state  at  8.92  eV,  with  the  next  E  state  at 
12.93  eV.  There  are  five  triplet  states  in  the  range  of  6  to  9  eV  for  the  D2(f  form, 
with  the  lowest  at  6.26.  The  triplet  states  start  at  7.07  eV  for  the  D3h  isomer,  with 
an  E'  state  at  7.60.  However,  despite  extensive  effort,  when  optimizing  the 
geometry  for  excited  states  either  with  CCSD  when  applicable,  or  EOM-CCSD 
analytical  gradient  techniques  [72]  to  determine  if  there  were  bound  excited  states, 
only  one  singlet  state  was  found.  Neither  the  energy  (shown  in  Table  2-1)  nor  the 
frequencies  (shown  in  Table  2-2)  make  a  persuasive  case  for  this  being  the  possible 
intermediate  state  in  the  1+1  experiment.  Taking  all  these  points  into 
consideration,  covalently  bound  energetic  tetraoxygen  molecules  do  not  appear 
likely  to  be  responsible  for  the  experimental  observations. 
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We  now  consider  electronically  excited  van  der  Waals  complexes.  There  are 
several  metastable  states  of  O2  that  may  form  long-lived  van  der  Waals  complexes 
of  the  correct  energy.  The  relevant  potential  energy  curves  for  O2  are  shown  in 
Figure  2-4.  Complexes  involving  lower  lying  metastable  singlet  states  are  known, 
but  have  neither  sufficient  energy  nor  plausible  ionization  paths  to  be  responsible 
for  the  experimental  observations.  The  Herzberg  states  of  02,  however,  do  appear 
at  the  correct  energy  to  form  complexes  that  could  give  rise  to  the  observed  spectra. 
One  of  these,  the  c  (1S~)  state,  is  shown  in  Figure  2-4  based  on  calculations  of 
Saxon  and  Liu  [74],  adapted  from  van  der  Zande  et  al  [73,75].  This  is  the  only 
species  that  may  possess  an  allowed  optical  transition  in  this  wavelength  region:  an 
excited  1  (1I19)  state  exists  very  near  the  energy  of  our  probe  transition  (i.e.,  about 
4  eV  above  this  c  state).  However,  earlier  calculations  indicated  that  this  1  (lUg) 
state  is  repulsive,  adiabatically  correlating  with  two  ground  state  oxygen  atoms.  In 
an  illuminating  series  of  experiments  [73, 75]  studying  atomic  fragments  after  charge 
transfer  from  cesium  atoms  to  O^",  van  der  Zande  et  al.  explored  the  nonadiabatic 
dynamics  and  coupling  among  several  of  the  curves  in  Figure  2-4  (and  with  triplet 
curves  not  shown).  Most  importantly,  they  showed  that  the  1  (lHlg)  and  2  ^Ug) 
curves  shown  in  Figure  2-4  do  not  interact  strongly,  and  may  be  viewed  in  a 
'diabatic'  picture,  as  shown;  the  1  (1119)  level  is  a  bound  or  quasibound  state.  The  1 
(^Ug)  <—  c  (XS~)  transition  of  02  is  thus  a  strongly  allowed  optical  transition  in  02 
with  nearly  diagonal  Franck-Condon  factors  in  the  region  of  280-330  nm.  It  is  not 
clear  that  this  important  implication  of  the  observations  of  van  der  Zande  et  al.  had 
been  recognized. 

However,  direct  ionization  of  this  1  (xns)  state  is  not  possible  in  this  wavelength 
region,  since  the  Franck-Condon  factors  connecting  it  with  the  ground  state  of  the 
ion  are  negligible.  In  fact,  the  experiments  of  van  der  Zande  and  coworkers  show  the 
importance  of  the  interactions  involving  the  1  (1n9)  valence  state  and  the  d  (XIIS) 
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Rydberg  state  that  was  initially  prepared  in  their  experiments.  This  shows  a  path 
to  ionization  in  this  wavelength  region  via  the  valence-Rydberg  interactions.  These 
considerations  yield  the  following  scenario  for  this  1+1  ionization  process  in  O4, 
indicated  by  the  heavy  arrows  in  Figure  2-4.  If  we  begin  with  a  van  der  Waals 
complex  between  O2  X  (3£~)  and  O2  c  (1S~),  a  fully  allowed  electronic  transition 
localized  on  the  c  state  molecule  takes  us  to  a  complex  involving  the  1  (1ng)  state. 

This  state  can  either  predissociate  to  give  oxygen  atoms  (and  an  O2  molecule), 
or  couple  to  the  d  (1ns)  Rydberg  state;  or  the  system  can  dissociate  to  two  O2 
molecules.  It  is  likely  that  all  of  these  occur,  no  doubt  with  a  strong  dependence  on 
the  initially  excited  vibrational  level.  If  the  Rydberg  complex  is  formed,  it  can  then 
ionize  easily  in  this  wavelength  region,  and  the  ionization  will  be  dominated  by 
Av  =  0  transitions  owing  to  the  diagonal  Franck-Condon  factors  between  the 
Rydberg  and  the  ion.  If  this  picture  is  accurate  for  O4,  it  is  surprising  that  no  O^"  is 
seen;  this  implies  some  significant  differences  for  the  ionization  dynamics  in  the 
complex  as  opposed  to  the  free  O2  molecule.  In  fact,  it  is  precisely  in  the  nature  of 
these  Rydberg-valence  interactions  that  we  can  expect  a  profound  impact  of  the 
formation  of  the  van  der  Waals  complex.  This  is  because  the  Rydberg  state  will  be 
greatly  stabilized  in  the  complex — nearly  to  the  extent  of  the  0.45  eV  bond  in 
C^-O^ .  The  valence  state  curves  will  be  little-perturbed  in  comparison.  The 
location  of  the  Rydberg  and  ion  curves  for  the  complex  are  shown  as  dashed  lines  in 
Figure  2-4.  This  provides  a  reasonable  explanation  for  the  absence  of  the  O^"  in 
these  experiments  despite  the  likelihood  that  the  number  density  of  free  02  c  (1E~) 
molecules  is  much  greater  than  those  involved  in  complexes.  The  fate  of  the  free  O2, 
upon  excitation  to  the  1  (^Tlg)  state,  is  either  predissociation  via  the  2  (1ns)  state, 
or  by  the  triplet  states  interacting  with  the  d  state. 

Many  of  the  experimental  results  can  be  satisfactorily  accounted  for  by 
invoking  this  complex.  The  rotational  spacing  in  the  long  wavelength  region,  about 
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3.5  cm-1,  is  very  near  4Be  for  the  c  (1S~)  state  (Be  =  0.9  cm-1).  This  would  be 
expected,  for  example,  for  a  T-shaped  complex  wherein  one  of  the  rotational 
constants  will  resemble  that  of  one  of  the  O2  molecules.  The  photoelectron  spectra, 
dominated  by  single  peaks,  arise  because  the  ionization  takes  place  from  a  complex 
involving  the  d  (1ns)  Rydberg  state  so  that  Av  =  0  transitions  dominate  as 
mentioned  above.  Finally,  the  absence  of  O^  is  readily  explained  by  the  very 
different  Rydberg- valence  interactions  in  the  complex  as  opposed  to  the  free  02. 
This  picture  also  accounts  for  some  unusual  spectra  reported  in  a  closely  related 
study  by  Helm  and  Walter  [76] .  Their  experiments  were  similar  to  the  studies  of 
van  der  Zande  et  al.,  but  used  charge  transfer  to  O4"  rather  than  Oj.  They  reported 
clearly  resolved  vibrational  structure  in  the  02  product  kinetic  energy  distributions 
after  charge  transfer  from  cesium,  which  they  reluctantly  ascribed  to  coincident 
formation  of  two  02  molecules  in  v=29,  a  rather  unlikely  process.  This  was 
necessary  to  account  for  the  vibrational  spacing  of  800  cm-1  observed  in  the  02 
kinetic  energy  release  distributions.  Our  alternative  interpretation  of  their  results 
suggests  simply  the  reverse  of  the  ionization  process  outlined  above:  electron 
transfer  from  cesium  populates  the  Rydberg  state  around  7.6  eV,  which  then 
couples  efficiently  to  the  metastable  02  X  (3E~)— 02  1  ^Ug)  complex.  We  suggest 
that  the  structure  in  the  kinetic  energy  release  distributions  of  Helm  and  Walter 
simply  reflects  the  vibrational  structure  in  the  metastable  state.  For  the  Herzberg 
states,  the  vibrational  frequencies  are  all  on  the  order  of  800  cm-1;  the  vibrational 
frequency  in  the  1  ^Ug)  state  is  likely  to  be  similar. 

It  is  important  to  note  that  although  these  spectra  are  not  associated  with 
covalently  bound,  energetic  04  species,  they  may  well  be  present  in  the  molecular 
beam. 


CHAPTER  3 
STABILIZATION  OF  THE  PSEUDO-BENZENE  N6  RING  WITH  OXYGENa 

3.1     Introduction 

One  of  the  goals  for  the  development  of  new  high-energy-density  molecules  has 

been  to  explore  the  prospects  for  making  novel  polynitrogen  compounds.  Since  the 

generation  of  N2  as  a  propulsion  or  explosion  product  is  highly  recommended  by  its 

very  strong  triple  bond, 

citehedmabstract,lauderdale,ferris,korkin,nreview, polynitrogen, perera  have  used 

predictive  quantum  chemical  techniques  to  investigate  the  structure,  stability, 

spectra,  and  decomposition  paths  of  experimentally  unknown  polynitrogen 

molecules,  as  have  others  [15,25,77-79].  (See  Reference  [80]  for  an  extensive  survey 

of  pure  nitrogen  species,  from  2  to  12  atoms,  their  cations,  anions  and  low-lying 

excited  states).  It  is  apparent  that  if -CH  groups  could  be  replaced  by  isovalent  -N, 

there  will  be  a  large  increase  in  the  amount  of  energy  stored  by  virtue  of  the 

standard  state  of  nitrogen,  N2,  and  also  in  the  repulsion  of  electron  lone  pairs  on 

adjacent  nitrogens.  However,  we  pay  a  price  for  this  endothermicity,  since  to  retain 

the  energy  for  use  as  a  fuel,  we  must  have  sufficiently  high  activation  barriers  to 

unimolecular  decomposition  as  well,  and  we  have  also  to  address  the  question  of 

whether  non-radiative  transitions  can  play  a  role  in  the  molecule's 

decomposition  [25].  We  have  considered  the  molecule  N4  in  a  tetrahedral 

arrangement  [21,24,28],  as  well  as  the  N8  analog  of  cubane  [24,81].  In  both  cases, 


a  This  chapter  consists  of  an  article  reproduced  with  permission  from  Kenneth  J. 
Wilson,  S.  Ajith  Perera,  Rodney  J.  Bartlett,  and  John  D.  Watts,  Journal  of 
Physical  Chemistry  A  volume  105,  7693-7699.  ©2001  American  Chemical  Society. 
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the  high  symmetry  indicates  that  any  decomposition  to  N2  would  be 
Woodward-Hoffman  forbidden,  suggesting  significant  barriers  to  decomposition. 
Numerical  investigation  of  the  barriers  shows  that  N4's  is  62  kcal  mol-1  [27,77] 
though  its  effective  barrier  is  close  to  30  (because  of  a  low-lying  triplet  state  [25,27]) 
and  N8's  is  19  [15,79].  In  Chapter  4  of  this  dissertation,  we  have  also  investigated 
the  other  azacubanes  where  1-7  of  the  -CH  units  have  been  replaced  by  -N,  and  the 
nitroazacubanes.  More  energy  per  molecule  can  be  gained  from  single  rather  than 
double-bonded  nitrogen  linkages,  recommending  N4  and  N8  if  they  could  be 
synthesized  and  stabilized,  but  at  the  cost  of  the  stability  offered  by  the 
double-bonded  -N=N-  structures.  In  addition,  to  impose  some  organization  of 
proposed  energetic  species,  we  here  consider  dimers,  trimers,  etc.  of  highly  energetic 
units  such  as  O2NCN  and  other  units  [82].  In  this  paper,  the  unit  is  N20. 

Beyond  N2  and  N3  ,  there  is  little  experimental  evidence  to  date  that  nitrogen 
can  form  stable  homonuclear  molecules.  Just  recently,  Christe  and  coworkers  have 
reported  the  synthesis  of  N^AsF^  [37]  while  Nj  [83],  N|  [84-86]  and  the  diazidyl 
N^"  complex  [87, 88]  have  been  detected  spectroscopically  as  short-lived  species. 
Notable  in  the  series  of  homoleptic  polynitrogen  systems  is  the  absence  of  the  N6 
ring.  The  N6  ring,  analogous  to  benzene,  is  not  a  minimum  on  the  potential  energy 
surface  as  a  planar,  hexagonal  ring.  Instead,  including  electron  correlation,  the  D6/j 
structure  is  a  second-order  saddle  point  with  the  closest  minimum,  a  non-planar 
boat- type  D2  structure  [89].  This  feature  might  be  attributed  to  the  large  lone-pair 
repulsion  in  the  ring.  However,  we  might  expect  to  reduce  some  of  this  repulsion  by 
forming  coordinate  covalent  bonds,  -N— >•().  In  fact,  such  bonds  are  found  in 
explosives  [90-93]  and  suggest  some  interesting  consequences  for  the  formation  of 
novel  polynitrogen  rings. 
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3.2     Computational  Methods 

All  the  ab  initio  calculations  in  this  report  were  performed  with  the  ACES  II 
program  system  [68]  while  the  density  functional  calculations  were  performed  with 
the  Q-CHEM  computer  program  [94].  The  well-known  hierarchy  of  ab  initio 
methods  ranging  from:  Hartree-Fock  self-consistent  field  (SCF),  second-order 
many-body  perturbation  theory  (MBPT(2)),  coupled-cluster  singles  and  doubles 
(CCSD)  to  coupled-cluster  singles  and  doubles  with  a  perturbative  inclusion  of 
triples  (CCSD(T))  were  employed  to  show  the  importance  of  higher  dynamic 
correlation.  Density  functional  theory  with  the  empirical  Becke  three-parameter 
Lee,  Yang  and  Parr  (B3LYP)  exchange  correlation  functional  was  also  used.  Both 
the  ab  initio  and  DFT  calculations  were  performed  in  the  6-31G*  one-particle  basis 
set  using  Cartesian  Gaussians  [69].  This  basis  corresponds  to  a  modest  15  functions 
per  heavy  atom,  but  is  largely  sufficient  for  the  whole  series  of  molecules  ranging  up 
to  nine  heavy  atoms,  studied  here.  To  obtain  better  energetics,  single  point  energies 
were  computed  with  the  aug-cc-pVDZ  basis  set  [95]  at  the  optimized  6-31G* 
geometries  (see  Table  3-1).  Total  atomic  charges  were  computed  with  the  natural 
bond  orbital  (NBO)  procedure  of  Weinhold  et  al  [96].  A  spin  unrestricted  wave 
function  was  used  to  describe  the  X  (3S~)  state  of  02  which  showed  little  evidence 
of  spin  contamination.  All  transition  states  were  confirmed  by  the  presence  of  only 
one  negative  eigenvalue  in  the  computed  Hessian  matrix.  All  core  electrons  were 
omitted  from  the  correlation  procedure,  however,  the  virtual  orbitals  corresponding 
to  the  core  electrons  where  included  in  the  correlated  calculation. 

3.3     Results  and  Discussion 

In  this  section,  we  begin  by  checking  the  reliability  of  our  methods  by 
calculations  on  the  N20  molecule.  Then,  we  explore  the  trends  in  structures  and 
energetics  for  its  dimers  and  trimers  plus  related  species  as  a  function  of  increasing 
the  number  of  coordinate-covalent  bonds.  The  first  series  is  based  on  a 
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Figure  3  1.  Computed  and  experimental  structure  for  N20. 

four-membered  nitrogen-ring  system  and  contains:  N40  and  N4O2.  In  the  second 
series,  we  consider  the  six-membered  nitrogen-ring  systems:  N60,  N6C"2  and  N603. 

The  computed  and  experimental  geometries  of  N20  are  shown  in  Figure  3-1. 
The  SCF  method  slightly  underestimates  both  bond  lengths,  while  all  other 
methods  predict  bond  lengths  longer  than  their  experimental  values.  The  computed 
CCSD(T)  N-N  bond  length  of  1.148  A  and  N-0  bond  distance  of  1.205  A  compare 
to  the  experimental  values  of  1.127  A  and  1.185  A,  respectively.  The  error  is  due  in 
large  part  to  the  use  of  the  small  6-31G*  one-particle  basis  set,  as  the  larger 
cc-pVTZ  basis  gives  1.132  A  and  1.189  A.  However  the  predicted  B3LYP  bond 
lengths  seem  slightly  better  in  the  6-31G*  small  basis.  The  computed  and 
experimental  harmonic  frequencies  are  presented  in  Table  3-2.  Here,  the  CCSD(T) 
computed  frequencies  are  slightly  better  than  the  B3LYP  computed  values.  The 
largest  difference  is  for  the  second  £+  normal  mode  which  CCSD(T)  underestimates 
by  6  cm-1  and  the  B3LYP  method  overestimates  by  88  cm-1. 

The  first  member  of  the  N4  series  is  the  N40  system  and  its  computed 
lowest-energy  structure  is  shown  in  Figure  3-2.  The  four-membered  nitrogen  ring  is 
not  a  minimum  for  this  molecule  which  instead  adopts  a  C3l,  structure  similar  to 
that  of  the  tetrahedral  isomer  of  N4.  One  of  the  N-N  bond  lengths  is  on  the  order  of 
a  single  bond  (1.45  A  [99])  with  its  CCSD(T)  predicted  value  of  1.458  A  and  the 
other  one  that  is  not  equivalent  by  symmetry  has  a  slightly  longer  value  of  1.536  A. 
Other  isomers  of  N40  have  been  considered  including  a  C2„  ring  structure  with 
oxygen  as  part  of  the  ring  [100].  However,  the  transition  state  for  the  unimolecular 
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Figure  3-2.  Computed  structure  for  N40. 
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Figure  3-3.  Computed  transition  state  structure  for  N40. 


Figure  3-4.  Lewis  structure  of  £rans-nitrosyl  azide. 


dissociation  of  this  C2v  molecule  was  located  which  indicated  a  barrier  between  1 
and  2  kcal  mol"1.  A  somewhat  more  stable  and  experimentally-observed  [101] 
isomer  is  the  toms-nitrosyl  azide  shown  in  Figure  3-4.  Prompted  by  its  Raman 
characterization,  Klapotke  and  Schulz  considered  two  dissociation  pathways  for  this 
molecule  [102].  The  first  pathway  involved  conversion  into  a  cis  isomer,  then  a  cyclic 
form,  followed  by  its  dissociation  into  N2  and  linear  N20.  The  highest  barrier  for 
this  process  computed  at  the  MP2/6-31+G*  level  was  6.7  kcal  mol-1.  A  transition 


24 

Table  3  3.  Activation  energies  in  kcal  mol-1  for  unimolecular  dissociation.  Values 
not  in  parentheses  are  from  the  aug-cc-pVDZ  basis  at  the  6-31G* 
optimized  geometries.  All  values  include  electronic  energy  differences  as 
well  as  zero-point  vibrational  energy  differences  and  thermodynamic 
corrections  for  298  K. 
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Figure  3-5.  Computed  structure  for  N402. 


state  for  the  direct  dissociation  into  N2  and  cyclic  N20  was  also  found  which 
indicated  a  barrier  of  24.2  kcal  mol-1. 

The  transition  state  structure  for  the  loss  of  NO  from  our  C3„  structure  is 
shown  in  Figure  3-3.  It  has  Cs  symmetry  with  the  major  difference  that  one  of  the 
N-N  bonds  is  broken.  It  is  worth  noting  that  the  differences  between  the  MBPT(2), 
CCSD,  CCSD(T)  and  B3LYP  methods  are  slightly  larger  than  those  for  the 
structure  shown  in  Figure  3-2,  indicating  that  more  sophisticated  treatments  of 
electron  correlation  are  needed  to  describe  transition  states.  The  barriers,  shown  in 
Table  3-3,  differ  at  most  by  7  kcal  mol"1  from  our  most  reliable  estimate,  the 
CCSD(T)  value  of  5.6  kcal  mol-1.  Such  a  small  value  suggests  that  it  may  be 
possible  to  observe  this  N40  isomer  only  as  a  short-lived  species,  rather  than  being 
suitable  for  preparation  and  handling  in  bulk  quantities. 
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Figure  3-6.  Computed  transition  state  structure  for  N402. 

The  addition  of  a  second  oxygen  atom  to  the  N4  ring  leads  to  the  completely 
planar  N402  structure  shown  in  Figure  3-5.  The  optimized  SCF  structure  belongs 
to  the  D2h  point  group  while  all  other  methods  predict  a  less  symmetric  C2/j 
structure.  One  N-N  bond  length  is  between  that  of  a  single  (1.45  A)  and  double 
(1.25  A)  bond  and  the  other,  not  equivalent  by  symmetry,  is  slightly  longer  than  a 
typical  single  bond.  Previously,  Manaa  and  Chabalowski  have  considered  another 
cyclic  isomer  of  N402  where  the  oxygen  atoms  were  members  of  the  ring  [103].  This 
structure  does  not  allow  for  oxygen  to  remove  as  much  charge  from  the  nitrogen 
atoms  and  hence  is  not  as  stable  kinetically.b    For  the  linear  N402  isomer,  there 
have  been  two  reports  of  its  synthesis  [104, 105]. 

The  addition  of  the  second  coordinate  covalent  bond  to  the  four-membered  ring 
significantly  increases  its  kinetic  stability.  The  transition  state  structure  for  the  loss 
of  N20  is  presented  in  Figure  3-6.  This  structure  has  Cs  symmetry  and  two 
considerably  longer  N-N  bonds.  The  large  differences  between  the  SCF  and  CCSD 


b  The  barrier  for  the  unimolecular  dissociation  of  the  N402  boat  isomer  is  11.0 
kcal  mol"1  at  the  CCSD(T)/6-3lG*  level.  Watts,  J.  D.;  Wilson,  K.  J.;  Bartlett,  R. 
J.  unpublished  work. 
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Figure  3-7.  Computed  structure  for  N60. 
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Figure  3-8.  Computed  transition  state  structure  for  N60. 

structures  demonstrate  the  importance  of  electron  correlation  in  describing  the 
transition  state.  Furthermore,  there  is  a  14  kcal  mol-1  difference  in  activation 
energies  with  the  SCF  computed  value  of  60.1  kcal  mol-1  and  the  CCSD  value  of 
45.6  kcal  mol-1.  Because  of  large  T2  amplitudes,  we  were  unable  to  converge  the 
CCSD(T)  method  on  this  transition  state. 

The  N60  molecule  is  the  first  member  of  the  second  series  with  a  six-membered 
ring.  Its  structure  is  nonplanar  with  Cs  symmetry  as  shown  in  Figure  3-7.  All  N-N 
bond  lengths  are  greater  than  that  of  a  double  bond  (1.25  A),  but  shorter  than  a 
single  bond  (1.45  A  [99]).  One  interesting  aspect  of  the  frequencies  and  intensities 
in  Table  3-6  is  the  substantial  difference  in  the  B3LYP  and  CCSD(T)  intensities. 
Particularly,  the  three  a'  modes  near  1000  cm-1  have  a  different  distribution  of 
intensities,  though  the  total  intensity  in  that  symmetry  is  comparable. 
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Figure  3-10.  Computed  transition  state  structure  for  N602- 


The  transition  state  for  the  unimolecular  dissociation  of  N60  is  shown  in  Figure 
3-8  which  corresponds  to  loss  of  N2.  Our  best  estimate  of  the  barrier  for  this  process 
is  1.2  kcal  mol"1  at  the  CCSD(T)/aug-cc-pVDZ//CCSD(T)/6-3lG*  level  of  theory. 
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Figure  3-11.  Computed  structure  for  N603. 


Addition  of  another  oxygen  atom  to  the  previous  structure  results  in  the  planar 
N602  system  shown  in  Figure  3-9.  This  molecule  has  C2,,  symmetry  and  all  N-N 
bond  lengths  are  between  those  of  a  double  and  single  bond.  Again,  the  IR 
intensities,  shown  in  Table  3-6  vary  greatly  between  CCSD(T)  and  B3LYP. 

The  transition  state  for  unimolecular  dissociation  is  shown  in  Figure  3-10  and 
corresponds  to  loss  of  N3.  The  barrier  for  this  process  is  13.8  kcal  mol-1  at  the 
CCSD(T)/aug-cc-pVDZ//CCSD(T)/6-31G*  level.  Although  this  value  is  somewhat 
larger  than  that  for  N60,  it  is  still  too  low  to  facilitate  handling  on  a  bulk  scale. 
Other  transitions  might  occur  to  products  involving  N^~  and  N30+,  the  latter 
isovalent  to  N4. 

Addition  of  another  oxygen  results  in  the  highly  symmetric  N603  structure 
shown  in  Figure  3-11,  the  trimer  of  N20.  This  molecule  has  D3h  symmetry  with  all 
N-N  bond  lengths  equal  to  1.363  A  at  the  CCSD(T)  level  of  theory.  This  distance  is 
almost  exactly  between  that  of  a  double  and  single  bond,  attesting  to  benzene-like 
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Figure  3-12.  Computed  transition  state  structure  for  N603. 


derealization  (though  no  resonance  stabilization  as  discussed  in  Section  3.5). 
Unlike  N60  and  N602,  N603  has  comparable  IR  intensity  patterns  between  B3LYP 
and  CCSD(T). 

The  transition  state  for  unimolecular  dissociation  is  shown  in  Figure  3-12.  Our 
best  estimate  at  the  CCSD(T)/6-31G*  level  for  the  barrier  of  this  process  is  62.4 
kcal  mol-1.  Being  the  trimer  of  N20,  there  could  be  an  alternative  dissociation  path 
to  3  N20's,  analogous  to  the  concerted  triple  dissociation  of  s-tetrazine  [106], 
however  all  attempts  to  locate  other  transition  states  along  different  dissociation 
pathways  were  unsuccessful.  Initial  indications  are  that  such  a  N603  species  should 
be  capable  of  being  synthesized. 

3.4     Atomic  Charges 

One  property  which  can  be  related  to  the  computed  activation  energies  is  the 
computed  atomic  charges  or  rather  the  arrangement  of  the  atomic  charges.  These 
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were  computed  with  the  natural  bond  orbital  (NBO)  formalism  [107]  within  the 
6-31G*  basis  and  are  presented  in  Figure  3-13.  In  all  structures,  oxygen  has  a  large 
negative  charge.  However,  in  the  stable  structures  (N4O2  and  N603),  the  nitrogen 
ring  is  composed  of  alternating  charges.  For  example  in  N4O2,  the  smallest  charge 
difference  on  adjacent  nitrogens  is  0.60  at  the  B3LYP  level  of  theory.  In  N60  and 
N602,  the  charge  on  adjacent  nitrogens  is  not  as  well  separated  with  the  smallest 
B3LYP  differences  being  0.07  and  0.13,  respectively.  Figure  3-13  shows  a  perfectly 
alternating  NBO  charge  on  the  six-membered  ring  of  N603  with  the  smallest  charge 
difference  being  0.56. 

3.5     Resonance  Stabilization? 

Another  interesting  concept  with  the  proposed  highly  energetic  molecules  is 
their  resonance  energy  and  whether  they  are  resonance  stabilized.  Early 
calculations  on  planar-hexagonal  N6  based  on  Shaik's  quantum-mechanical 
resonance  energy  [108]  and  Dewar's  7r-resonance  energy  [109]  indicated  that 
hexazine  was  even  more  aromatic  than  benzene.  However,  Glukhovtsev  and 
Schleyer  later  reported  homodesmic  and  hyperhomodesmic  reactions  for  hexazine 
which  showed  a  destabilizing  resonance  energy  of  17.6  kcal  mol-1  and  10.4  kcal 
mol-1,  respectively  [110,111]. 

Gimarc  and  Zhao  have  offered  an  explanation  of  nitrogen's  destabilizing 
resonance  energy  as  opposed  to  that  of  carbon's,  which  is  based  on  average  bond 
energies  [112,113].  In  this  approximation,  the  total  energy  of  a  molecule  is  the  sum 
of  its  bond  energies.  Since  the  average  carbon-carbon  single  bond  energy  is  83  kcal 
mol-1  which  is  more  than  half  of  the  carbon-carbon  double  bond  energy  of  144  kcal 
mol-1,  carbon  structures  with  pairs  of  single  bonds  rather  than  double  bonds  are 
lower  in  energy.  For  nitrogen,  the  situation  is  reversed.  The  average 
nitrogen-nitrogen  single  bond  energy  is  43  kcal  mol"1,  which  is  less  than  half  of  the 
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nitrogen-nitrogen  double  bond  energy  of  100  kcal  mol-1.  Hence,  nitrogen  prefers 
double  bonds  rather  than  pairs  of  single  bonds. 

To  estimate  the  amount  of  resonance  energy  in  the  novel  N4O2  and  N603 
molecules,  we  have  used  the  following  isodesmic  reactions,  in  analogy  to  N6  +  3 
NH3  — *  3  H2N-N=NH. 


h(  .n        +2NH3  ►         2  N=NH 

\v  / 

N  HjN 

I 

O 


O 
.3NH3  ►         3         \=NH 


O 


^•N\  •*< 


Figure  3-14.  Isodesmic  reactions  for  N402  and  N603. 

Based  on  the  energies  of  the  isodesmic  reactions  shown  in  Figure  3-14,  both  ring 
structures  have  a  destabilizing  resonance  energy.  For  N4O2,  the 
CCSD(T)/aug-cc-pVDZ//CCSD(T)/6-31G*  value  is  27.6  kcal  mol"1  and  for  N603, 
13.9  kcal  mol-1.  Other  quantum  mechanical  methods  predict  values  which  agree 
closely  with  CCSD(T)  and  they  are  presented  in  Table  3-4.  We  also  present  the 
resonance  energy  for  N6  in  the  D6h  conformation.  In  fact,  the  value  of  -17.6  is  very 
close  to  our  MBPT(2)/6-3lG*  value  for  N603  of -18.6  kcal  mol  -1.  For  N402,  the 
resonance  energy  is  slightly  more  destabilizing,  numerically  being  -36.4  kcal  mol-1. 
3.6     Enthalpy  of  Formation  and  Specific  Impulse 
One  useful  figure  of  merit  for  potential  fuels  and  explosives  is  the  material's 
enthalpy  of  formation,  AHy,  or  energy  relative  to  the  elements  in  their  standard 
states.  It  has  been  shown  for  molecules  that  do  not  contain  fluorine,  the  enthalpy  of 
formation  largely  parallels  the  heat  of  combustion  [115].  In  Table  3-5,  we  present 
AH/'s  for  the  ring  molecules  considered  in  this  work  and  the  N20  test  system.  For 
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Table  3-5.  Enthalpy  of  formation  (AH/)  and  specific  impulse  (lsp)  values  for  the 
ring  species  considered  in  this  work.  AH/  is  in  units  of  kcal  mol-1  and 
Isp  is  in  seconds.  Values  not  in  parentheses  are  from  the  aug-cc-pVDZ 
basis  at  the  6-31G*  optimized  geometries.  All  values  include  electronic 
energy  differences  as  well  as  zero-point  vibrational  energy  differences  and 
thermodynamic  corrections  for  298  K. 


AH7 

'■sp 

B3LYP 

experiment  [114] 

B3LYP 

experiment0 

N20 

16.7  (17.3) 

19.61 

163.2  (166.3) 

176.9 

N40 

205.2  (206.4) 

447.3  (448.6) 

N402 

125.4  (125.0) 

316.3  (315.8) 

N60 

169.1  (169.3) 

344.5  (344.7) 

N602 

160.4  (160.2) 

311.5  (311.4) 

N603 

155.6  (154.7) 

287.7  (286.8) 

a.)  Calculated  using  the  experimental  AH/  and  Equation  3-1. 

N20,  the  AH/  computed  from  B3LYP/aug-cc-pVDZ//B3LYP/6-31G*  is  17.3  kcal 
mol"1  which  is  in  excellent  agreement  with  the  experimental  value  of  19.61  kcal 
mol-1  [114].  For  propellants,  the  molecular  weight  is  also  important  and  a 
material's  potential  is  best  measured  by  its  specific  impulse,  lsp.  The  specific 
impulse  in  units  of  seconds  can  be  approximated  with  the  equation  [116]: 


-   A/// (kcal  mol"1)  .       . 

/sp(seconds   =  265  J  — — f -J-  3-1 

''  MW (grams  mol  x) 

where  M W  is  the  molecular  weight  in  grams  per  mol.  The  L,p's  for  the  ring  species 
considered  in  this  work  are  presented  in  Table  3-5.  The  prospective  HEDMs  that 
are  stable  with  respect  to  unimolecular  dissociation,  N402  and  N603,  have  Isp's  of 
315.8  and  286.8  seconds,  respectively.  Both  of  these  offer  an  improvement  to  the 
224  seconds  Isp  for  hydrazine,  the  most  frequently  used  monopropellant.  A  survey  of 
possible  generalizations  of  the  basic  hydrazine  molecular  structure  has  been 
considered  elsewhere  [117]. 
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Table  3-6.  Computed  normal  harmonic  vibrational  frequencies  and  infrared 

intensities  in  parentheses  for  various  nitrogen  ring  structures  within  the 
6-31G*  basis  except  for  the  CCSD(T)  frequencies  of  the  N603  isomer 
which  were  computed  with  the  larger  cc-pVDZ  basis  set.  Frequencies  are 
in  units  of  cm-1,  intensities  are  in  km  mol-1. 


SCF       MBPT(2) 


CCSD 


CCSD(T) 


B3LYP 


N40 


■3r 


e 
ai 

c 
ai 

ai 


522 

841 

834 

1173 

1378 

1726 


382 
558 
710 
762 
912 
1703 


465 

659 

772 

932 

1156 

1603 


430(2.4) 
584(6.7) 
730(1.1) 
845(15.2) 
1059(1.8) 
1579(344.7) 


453(2.9) 
604(10.0) 
778(2.3) 
884(15.7) 
1149(1.1) 
1657(347.3) 


N40  TS 


732i 


683i 


5991  (38.8)        612i  (49.1) 


N4O2 


'2h 


a9 

K 


251 

443 

575 

782 

805 

894 

931 

969 

1260 

1450 

1931 

2217 


224 

372 

558 

656 

690 

908 

737 

866 

1183 

1346 

1795 

1879 


228 

372 

519 

666 

714 

898 

772 

853 

1130 

1356 

1726 

1935 


222(2.2) 
328(4.8) 

486 

622 

702 

870 
704(14.6) 
789(93.6) 

1147 
1292(57.8) 
1673(617.8) 

1865 


228(2.1) 
343(4.2) 

498 

664 

725 

876 
743(19.5) 
813(102.0) 

1133 
1315(74.0) 
1706(758.0) 

1904 


N4O2  TS 


1063i 


14171  1070i  (73.4) 


801i  (23.5) 


N60 


a 

a' 

a' 

a" 

a' 

a" 

a' 

a' 

a' 

a" 

a" 

a' 

a' 


130 

451 

494 

690 

693 

675 

924 

793 

965 

1167 

1294 

1493 

1589 

1715 

1835 


153 
322 
490 
516 
669 
673 
752 
827 
857 
1055 
1138 
1166 
1262 
1299 
1717 


114 
260 
358 
549 
665 
574 
789 
710 
859 
1091 
1115 
1266 
1340 
1440 
1672 


136(0.5) 

313(2.9) 

425(19.1) 

514(17.3) 

651(0.7) 

569(32.2) 

746(4.5) 

692(5.5) 

821(19.1) 

1049(12.3) 

1073(23.7) 

1173(0.0) 

1257(11.6) 

1340(31.2) 

1609(201.6) 


148(0.5) 

318(1.0) 

493(5.1) 

532(19.2) 

674(1.0) 

573(24.7) 

774(2.6) 

711(7.3) 

835(30.4) 

1066(30.1) 

1091(2.0) 

1172(1.4) 

1284(23.6) 

1390(38.4) 

1648(269.0) 


N60  TS         Ci 


882i 


649i 


511i 


416i  (51.5)        448i  (52.9) 


3G 


Table  3-6 — continued. 


SCF        MBPT(2)        CCSD         CCSD(T) 


B3LYP 


N60 


6^2 


-'2r 


a2 
bi 
a, 
bi 
b2 
ai 
b2 
b2 
a2 
bi 
ai 
ai 
1)2 

Hi 

ai 

1)2 

b2 
ai 


168 
194 
518 
569 
663 
672 
703 
811 


877 

861 

1092 

1249 

1443 

1484 

1664 

1812 

1909 


98 

145 

439 

462 

541 

591 

605 

706 

714 

722 

768 

1041 

1085 

1211 

1294 

1413 

1653 

1681 


N602  TS         C, 


a 


474i 


541i 


N60 


6^3 


D 


3A 


e 

a2 

e' 
e' 

Bi' 

a2' 
a2" 

e" 
ai' 
a2' 

e' 

e' 

e' 
ai' 


166 

236 

512 

633 

791 

822 

830 

840 

1038 

1134 

1252 

1510 

1788 

1930 


131 

207 

442 

575 

674 

628 

670 

689 

988 

957 

1092 

1364 

1645 

1687 


113 

149 

458 

482 

489 

600 

612 

643 

747 

751 

754 

1010 

1088 

1259 

1295 

1451 

1628 

1702 


131 

206 

456 

572 

697 

648 

704 

721 

956 

916 

1102 

1346 

1610 

1720 


98 
141(0.4) 
438(4.1) 
453(16.0) 
486(17.5) 
578(1.7) 
595(1.9) 
616(48.8) 

709 

714(14.3) 

713(2.2) 

992(11.8) 

1036(30.1) 

1193(1.6) 

1226(22.3) 

1367(133.4) 

1565(447.0) 

1633(65.6) 


127 

201(3.3) 

442(7.6) 

552(12.9) 

662 

602 

669(33.7) 

701 

938 

875 

1052(41.3) 

1283(175.0) 

1559(469.5) 

1660 


90 
146(0.4) 
445(2.9) 
458(17.9) 
536(13.4) 
583(2.8) 
606(33.3) 
628(13.4) 

738 

745(18.2) 

706(4.2) 

986(7.9) 

1060(36.3) 

1237(14.7) 

1200(22.0) 

1409(175.8) 

1596(522.6) 

1672(84.7) 


448i  446i  (94.9)  319i  (49.0) 


128 

202(4.5) 

444(5.7) 

562(14.7) 

679 

690 

695(38.6) 

712 

943 

915 

1064(34.7) 

1271(218.4) 

1572(508.9) 

1682 


N603  TS 


638i 


364i 


604i  4201  (146.8)        293i  (114.2) 


3.7     Conclusions 
In  this  work,  we  have  quantitatively  shown  how  four  and  six-membered 
nitrogen  rings  can  be  stabilized  by  coordinate  covalent  bonds  to  oxygen.  Other 
potentially  interesting  coordinate  covalent  structures  would  include  those  to  BH3. 
Our  analysis  is  based  on  locating  the  lowest  energy  transition  state  for  unimolecular 
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dissociation.  Although  we  believe  we  have  done  this  effectively,  we  realize  that  with 
large  molecules,  and  hence  high-dimensional  potential  energy  surfaces,  other 
decomposition  routes  may  exist.  Similarly-transformed  equation-of-motion 
(STEOM)  calculations  [71]  of  the  vertical  excitation  energies  at  the  singlet 
optimized  geometries  of  N4O2  and  N603  show  the  lowest  triplet  states  to  be  at  60.4 
and  56.1  kcal  mol-1,  respectively.  Since  both  of  the  triplets  are  high  in  energy, 
singlet-triplet  crossings  should  not  significantly  lower  the  activation  energies.  In 
light  of  the  stability  characteristics  of  the  N20  dimer  and  trimer,  the  tetramer  of 
N20  might  be  expected  to  have  similar  properties.  Our  preliminary  investigations 
show  that  Ng04  has  a  S4-type  structure.  Although  the  smallest  charge  difference  on 
adjacent  nitrogens  of  the  tetramer  is  0.54  at  the  B3LYP  level  of  theory,  the  structure 
is  highly  nonplanar.  We  have  not  investigated  its  transition  states  for  unimolecular 
dissociation,  and  therefore  cannot  comment  on  its  kinetic  stability.  All  told,  N402 
and  N603  are  highly  energetic  molecules  which  appear  to  be  stable.  They  are  indeed 
worth  attempts  to  synthesize.  To  facilitate  their  identification,  we  present  harmonic 
infrared  vibrational  frequencies  and  their  associated  intensities  in  Table  3-6. 


CHAPTER  4 
HEATS  OF  FORMATION  FOR  THE  AZACUBANES  AND 
NITRO-SUBSTITUTED  AZACUBANES 

4.1     Introduction 
The  cubane  system  was  first  synthesized  over  35  years  ago  by  Eaton  and 
Cole  [118].  In  light  of  cubane's  immense  strain  energy  (166  kcal  mol"1)  and  large 
positive  heat  of  formation  (148.7  kcal  mol"1),  cubane's  kinetic  stability  up  to  230°C 
is  quite  unique  [119,120].  This  is  in  large  part  due  to  its  highly  symmetric  structure 
which  makes  many  of  the  dissociation  pathways  Woodward-Hoffman  forbidden. 
Furthermore,  appreciable  geometry  changes  are  only  possible  if  two  C-C  bonds  are 
broken  simultaneously.  Another  feature  of  cubane's  caged  structure  is  its  high 
density  of  1.29  g  cm"3  [120]. 

Based  on  their  untypical  structure  and  properties,  cubane  and  its  derivatives 
have  emerged  as  outstanding  candidates  for  high-energy  density  materials  [21,24]. 
In  fact,  many  of  the  nitrosubstituted  cubanes  have  been  prepared  including 
octanitrocubane  [121-124].  Since  nitrogen  is  isoelectronic  with  -CH,  some  obvious 
derivatives  are  the  azacubanes  (CnN8_„H„)  and  nitroazacubanes  (CnN8_n(N02)„). 
Chemical  intuition  based  on  the  standard  state  of  nitrogen  suggests  that  the 
azacubanes  will  be  higher  in  energy  than  the  cubane  parent  molecule.  However, 
intuition  does  not  provide  a  means  for  estimating  the  magnitude  of  this  energy 
difference.  Since  experimental  data  on  the  thermochemical  properties  of  azacubanes 
is  nonexistent  (although  one  potential  precursor  has  been  made  [125])  a  theoretical 
investigation  of  these  properties  is  warranted. 

Estimating  the  relative  stabilities  of  various  azacubanes  having  a  fixed  number 
of  nitrogen  atoms  (n)  does  not  present  a  difficult  theoretical  problem  since  the  gross 
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structural  features  of  the  possible  isomers  are  similar,  resulting  in  an  approximate 
cancellation  of  errors  in  the  electronic  structure  calculations.  However,  the  energy 
difference  between  different  classes  of  azacubanes  (e.g.  between  the  n=6  and  n=4 
isomers  of  C„N8_„Hn)  is  considerably  more  difficult  to  predict.  Indeed,  the  most 
significant  quantity  to  those  interested  in  the  potential  use  of  azacubanes  as  fuels  is 
the  molar  heat  of  formation  (AH/)  defined  by  the  enthalpy  change  in  the  reaction: 

nH2  +  -^N2  +  nC   — *  CnN8-nHn  (4-1) 

and  a  similar  equation  exists  for  the  nitroazacubanes.  Direct  evaluation  of  AHf  by 
calculation  of  electronic  energies  for  all  species  involved  together  with  corrections 
for  zero-point  vibrational  energies  and  temperature  effects  (with  the  latter 
contributions  hereafter  referred  to  as  "thermodynamic  corrections" )  is  a  notoriously 
difficult  theoretical  problem  due  to  the  different  bonding  situations  present  in  the 
reactant  and  product  sides  of  the  chemical  equation.  For  example,  it  has  been 
demonstrated  that  calculation  of  the  heat  of  formation  of  gaseous  ammonia  by 
direct  calculation  of  NHZ,  JV2,  and  H2  produces  results  that  can  oscillate  wildly 
with  the  level  of  theory  and  choice  of  basis  set. 

However,  a  procedure  recommended  by  Pople  and  collaborators  over  three 
decades  ago  based  on  the  use  of  isodesmic  reactions  [126, 127]  provides  a  convenient 
means  for  determining  heats  of  formation.  An  isodesmic  reaction  is  defined  as  one 
in  which  the  types  of  all  chemical  bonds  are  conserved  in  the  course  of  the  reaction. 
A  trivial  example  of  an  isodesmic  reaction  is  the  conversion  of  the  n=6  azacubane 
with  nitrogens  separated  by  the  body  diagonal  to  that  in  which  the  nitrogens  are 
separated  by  a  face  diagonal.  In  both  structures,  there  are  6   C    C  single  bonds,  6 
C-  N  single  bonds,  and  4   C-  H  bonds.  Interest  in  isodesmic  reactions  stems  from 
the  fact  that  the  shortcomings  of  approximate  wavefunctions  and  properties 
calculated  from  them  are  intrinsically  related  to  the  chemical  environment  of  the 
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atoms  making  up  a  molecule.  As  a  result  of  the  structural  similarities  present  on 
both  sides  of  the  equation,  errors  in  absolute  energies  calculated  for  the  species 
involved  in  an  isodesmic  reaction  tend  to  cancel  when  the  overall  energy  change  for 
the  reaction  is  calculated.  Ideally,  one  chooses  an  isodesmic  reaction  where  AH;  for 
all  species  except  the  one  of  interest  have  been  accurately  established 
experimentally.  It  is  then  a  straightforward  matter  to  calculate  the  heat  of 
formation  of  the  unknown  species  from  a  combination  of  the  experimental  results 
and  the  theoretical  enthalpy  change  for  the  isodesmic  process.  Such  an  approach 
has  been  used  with  considerable  success  in  the  past  to  calculate  heats  of  formation 
for  a  wide  variety  of  molecules. 

Nevertheless,  the  concept  of  chemical  environment  tacitly  assumed  in  the 
preceding  paragraph  is  oversimplified.  For  example,  one  might  use  the  isodesmic 
reaction: 

I6CJ/4  +  CsHz  — >  12C2H6  (4-2) 

together  with  experimental  values  for  ethane  and  methane  to  estimate  the  heat  of 
formation  of  cubane  via  the  relation: 

AHf(C8Hs)  =  12AHf(C2He)  -  16AHf{CH4)  -  HR  (4-3) 

where  HR  is  the  calculated  enthalpy  change  for  the  isodesmic  reaction.  Note  that 
while  Equation  4-2  conforms  to  our  definition  of  an  isodesmic  reaction  (12   C-  C 
bonds  and  72   C-  H  bonds  on  both  sides  of  the  equation),  it  is  clear  that  the   C-  C 
bonds  in  ethane  are  different  from  those  in  cubane,  which  is  a  highly  strained 
system.  Therefore,  one  should  expect  the  approximate  cancellation  of  errors  in 
isodesmic  reactions  of  this  type  to  be  less  satisfactory  than  for  processes  such  as: 

C3H6  +  CH4  — ►  C2He  +  C2H4  (4-4) 
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Table  4-1.  Naming  conventions  for  the  azacubanes  and  nitroazacubanes  investigated 
in  this  study.  The  atomic  positions  are  shown  in  Figure  4-1.  For  the 
azacubanes,  the  substituent  of  the  carbon  atom  is  hydrogen  (X=H),  and 
for  the  nitroazacubanes  the  substituent  is  a  nitro  group  (X=N02). 


Abbreviated  Name 

Formula 

1 

2 

3 

4 

5 

6 

7 

8 

(CX)8 

C 

C 

C 

C 

C 

C 

C 

C 

1-aza 

(CX)7N 

N 

c 

C 

C 

c 

C 

c 

c 

1,3-diaza  (trans) 

(CX)6N2 

N 

c 

N 

c 

c 

c 

c 

c 

1,8-diaza  (opp) 

(CX)6N2 

N 

c 

c 

c 

c 

c 

c 

N 

1,2-diaza  (cis) 

(CX)6N2 

N 

N 

c 

c 

c 

c 

c 

c 

1,3,5-triaza 

(CX)5N3 

N 

c 

N 

c 

N 

c 

c 

c 

1,2,5-triaza 

(CX)5N3 

N 

N 

c 

c 

N 

c 

c 

c 

1,2,3-triaza 

(CX)5N3 

N 

N 

N 

c 

c 

c 

c 

c 

1,3,5,7-tetraaza 

(CX)4N4 

N 

c 

N 

c 

N 

c 

N 

c 

1,2,3,5-tetraaza 

(CX)4N4 

N 

N 

N 

c 

N 

c 

c 

c 

1,2,5,8-tetraaza 

(CX)4N4 

N 

N 

C 

c 

N 

c 

c 

N 

1,2,3,7-tetraaza 

(CX)4N4 

N 

N 

N 

c 

C 

c 

N 

c 

1,2,3,6-tetraaza 

(CX)4N4 

N 

N 

N 

c 

C 

N 

c 

c 

1,2,3,4-tetraaza 

(CX)4N4 

N 

N 

N 

N 

c 

c 

c 

c 

1,2,3,5,7-pentaaza 

(CX)3N5 

N 

N 

N 

c 

N 

c 

N 

c 

1,2,3,5,6-pentaaza 

(CX)3N5 

N 

N 

N 

c 

N 

N 

c 

c 

1,2,3,4,5-pentaaza 

(CX)3N5 

N 

N 

N 

N 

N 

c 

c 

c 

1,2,3,4,5,7-hexaaza  (trans) 

(CX)2N6 

N 

N 

N 

N 

N 

c 

N 

c 

1,2,3,5,6,8-hexaaza  (opp) 

(CX)2N6 

N 

N 

N 

c 

N 

N 

c 

N 

1,2,3,4,5,6-hexaaza  (cis) 

(CX)2N6 

N 

N 

N 

N 

N 

N 

c 

c 

1,2,3,4,5,6,7-septaaza 

(CX)N7 

N 

N 

N 

N 

N 

N 

N 

c 

N8 

N 

N 

N 

N 

N 

N 

N 

N 

which  could  be  used  to  determine  the  heat  of  formation  of  propylene  from  known 
values  for  methane,  ethane,  and  ethylene.  A  large  volume  of  experience  has 
demonstrated  that  enthalpy  changes  for  isodesmic  reactions  of  the  latter  variety  can 
be  calculated  at  even  low  levels  of  theory  such  as  the  self-consistent  field  (SCF) 
approximation  with  minimal  or  split-valence  basis  sets.  However,  for  reactions  such 
as  Equation  4-2  in  which  the  conserved  "bond  types"  are  less  similar  chemically, 
higher  levels  of  theory  should  be  used. 

In  this  chapter,  we  apply  isodesmic  reaction  strategies  to  study  the  absolute 
heats  of  formation  for  all  isomers  of  the  n=8  to  n=0  azacubane  and  nitroazacubane 
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systems.  Due  to  the  strained  nature  of  these  systems  and  their  unusual  bonding 
environments,  the  considerations  discussed  above  suggest  that  calculations  based  on 
a  low-level  ab  initio  approach  may  not  be  satisfactory.  Hence,  we  have  investigated 
the  sensitivity  of  the  predicted  heats  of  formation  by  performing  calculations  at 
levels  of  sophistication  ranging  from  the  simple  SCF  model  to  coupled-cluster  (CC) 
treatments  that  include  effects  of  triple  excitations.  The  purpose  of  the  present 
study  is  threefold.  In  addition  to  providing  what  we  believe  to  be  accurate 
predictions  for  the  formation  enthalpies  of  the  azacubanes  and  the  nitroazacubanes 
beyond  what  is  in  the  current  literature  [128-131],  the  systematic  study  of  the 
dependence  of  AHf  as  a  function  of  the  correlation  treatment  should  provide  some 
guidelines  for  investigations  of  related  molecules.  Finally,  we  compare  the  results 
obtained  here  with  those  predicted  by  semiempirical  molecular  orbital  approaches 
and  heats  of  formation  computed  directly  from  the  standard  states. 
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Figure  4-1.  Numbering  scheme  used  in  this  work. 


Table  4-2.  Model  compounds  for  isodesmic  reactions  and  their  experimental  heats 
of  formation  in  kcal  mol^1  [114]. 


molecule 

ch7~ 

CH3NH2 
CH3NO2 
C2H6 
NH3 

N(CH3)3 
N2H4 


AH7 
-17.9 
-5.4 
-19.3 
-20.0 
-11.0 
-5.7 
22.8 
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4.2     Methods 

The  first  step  in  calculating  the  heat  of  formation  is  to  determine  an  isodesmic 
reaction  for  each  molecule.  A  small  number  of  simple  molecules  were  chosen  as 
sources  of  different  types  of  bonds.  For  example,  for  a   C-  C  single  bond,  ethane 
was  used.  For  a  N-  C  single  bond,  methyl  amine  was  used.  All  of  the  model 
compounds  are  listed  in  Table  4-2  along  with  their  experimental  heats  of  formation. 
The  sole  restriction  is  that  only  molecules  with  experimentally  known  heats  of 
formation  in  the  ideal  gas  state  may  be  used.  This  set  of  molecules  was  used  to 
form  balanced  chemical  reactions  where  both  the  individual  atoms  and  the  number 
of  each  type  of  bond  were  balanced.  Isodesmic  reactions  for  the  azacubanes  and 
nitroazacubanes  are  given  in  Tables  4-3  and  4-6,  respectively.  One  consequence  of 
the  large  stoichiometric  coefficients  is  that  small  errors  in  the  computed  energy  for 
each  model  compound  are  multiplied  by  large  numbers. 

In  this  study,  CC  calculations  were  done  with  the  ACES  II  package  [68].  The 
semi-empirical  and  DFT  calculations  were  performed  with  the  GAUSSIAN  94 
package  [132].  The  calculations  were  done  using  the  Dunning  double  (  (DZ)  basis 
set  [133,134]  with  polarization  functions  from  correlated  calculations  [135].  The 
geometries  for  all  molecules  were  optimized  at  the  SCF  level  and  are  available  upon 
request  to  interested  parties.  Using  these  geometries,  the  energies  were  calculated  at 
the  AMI  [136],  MIND03  [137],  PM3  [138],  BLYP,  B3LYP,  SCF,  MBPT(2),  CCSD, 
and  CCSD(T)  [139]  levels.  Pure  spherical  harmonics  (i.e.,  5  d-type  functions)  were 
used  and  all  core  electrons  were  omitted  from  the  correlation  procedure.  We  did  not 
perform  CCSD  and  CCSD(T)  calculations  for  the  series  of  nitroazacubanes. 

Frequency  calculations  were  also  done  at  the  SCF  level,  and  the  frequencies 
obtained  were  used  to  verify  the  existence  of  a  structure  with  no  imaginary 
frequencies.  Moreover,  the  computed  frequencies  were  used  to  calculate  the 
zero-point  vibrational  energy  as  well  as  the  thermodynamic  corrections  for 
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Table  4-3.  Reaction  coefficients  for  the  isodesmic  reaction:  azacubane  +  a  CH4  4-  /3 
NH3  -*  7  N2H4  +  5  CH3NH2  +  e  C2H6 


azacubane 

a 

P 

7 

S 

e 

(CH)8 

0 

16 

0 

0 

12 

1-aza 

2 

14 

0 

3 

9 

1,3-diaza  (trans) 

4 

12 

0 

6 

6 

1,8-diaza  (opp) 

4 

12 

0 

6 

6 

1,2-diaza  (cis) 

4 

12 

1 

4 

7 

1,3,5-triaza 

6 

10 

0 

9 

3 

1,2,5-triaza 

6 

10 

1 

7 

4 

1,2,3-triaza 

6 

10 

2 

5 

5 

1,3,5,7-tetraaza 

8 

8 

0 

12 

0 

1,2,3,5-tetraaza 

8 

8 

2 

8 

2 

1,2,5,8-tetraaza 

8 

8 

2 

8 

2 

1,2,3,7-tetraaza 

8 

8 

3 

6 

3 

1,2,3,6-tetraaza 

8 

8 

3 

6 

3 

1,2,3,4-tetraaza 

8 

8 

4 

4 

4 

1,2,3,5,7-pentaaza 

10 

6 

3 

9 

0 

1,2,3,5,6-pentaaza 

10 

6 

4 

7 

1 

1,2,3,4,5-pentaaza 

10 

6 

5 

5 

2 

1,2,3,4,5,7-hexaaza 

(trans) 

12 

4 

(i 

6 

0 

1,2,3,5,6,8-hexaaza 

(opp) 

12 

4 

6 

6 

0 

1,2,3,4,5,6-hexaaza 

(cis) 

12 

1 

7 

4 

1 

1,2,3,4,5,6,7-heptaaza 

11 

2 

9 

3 

0 

N8 

16 

0 

12 

0 

0 

finite-temperature.  Given  this  data,  the  heats  of  formation  for  each  of  the 
azacubanes  and  nitroazacubanes  were  determined. 

Heats  of  formation  were  also  computed  directly  from  the  standard  states  of  the 
elements.  Again,  this  is  an  extremely  demanding  test  for  a  theoretical  method  as 
nearly  every  bond  in  the  molecule  is  broken.  To  correct  for  the  standard  state  of 
carbon,  we  used  the  experimental  heat  of  formation  of  169.9  kcal  mol-1  for  the  3P 
state  [114]. 

4.3     Results  and  Discussion 

The  heats  of  formation  for  the  series  of  azacubanes  are  presented  in  Table  4-4. 
The  series  is  arranged  in  terms  of  increasing  nitrogen  content  and  for  the 
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azacubanes  with  the  same  chemical  formula,  increasing  N-  N  bonds.  For  cubane, 
there  are  two  experimental  values  of  its  heat  of  formation.  In  1966,  a  value  of  148.7 
kcal  mol-1  was  obtained  directly  from  cubane  [140].  A  more  recent  value  of  159  kcal 
mol-1  was  estimated  indirectly  from  the  heat  of  combustion  of 
l,4-bis(methoxycarbonyl)cubane  [141].  Our  most  reliable  CCSD(T)  estimate  of 
150.4  kcal  mol-1  is  more  consistent  with  the  older,  direct  measurement.  Trends  in 
the  AHf  values  for  different  methods  are  better  presented  graphically  in  Figure  4-2. 
Here,  the  AMI  and  MIND03  methods  show  large  variations  from  all  of  the  other 
methods.  However,  PM3,  the  other  semiempirical  method,  is  closer  to  the  ab  initio 
methods.  Figure  4-3  offers  a  magnified  view,  where  AH/  differences  with  respect  to 
CCSD(T)  are  plotted.  The  BLYP  and  B3LYP  values  show  the  same  trend  with 
large  variations  along  the  series.  The  MBPT(2)  and  CCSD  results  more  closely 
follow  the  CCSD(T)  values  for  the  range  of  molecules  studied. 

Another  important  trend  along  the  series  of  azacubanes  is  the  energy  change  in 
the  isodesmic  reactions.  The  isodesmic  reaction  energies  for  the  various  methods  are 
presented  in  Table  4-5  and  graphically  in  Figure  4-4.  In  the  bond  energy  additivity 
model  where  the  energy  of  a  molecule  is  approximated  by  the  sum  of  its  bond 
energies,  these  values  would  be  zero  as  they  represent  strain,  resonance 
stabilization,  and  other  effects.  For  nearly  all  of  the  azacubanes,  -AEiSO(iesmic  is 
greater  than  zero  suggesting  that  they  are  less  stable  than  their  molecular 
fragments.  However,  -AEisodesmic  does  decrease  along  the  series  with  increasing 
nitrogen  content.  Beginning  with  cubane,  -AEisodesmic  is  104.6  kcal  mol-1  from 
CCSD(T)  and  decreases  to  -2.8  kcal  mol-1  for  N8.  This  trend  is  consistent  with  an 
earlier  investigation  which  showed  that  introduction  of  nitrogen  into  the  cube 
stabilizes  it  by  a-conjugation  of  the  lone  pairs  [142].  Also,  since  nitrogen  prefers 
angles  slightly  less  than  109.5  degrees  as  opposed  to  carbon,  more  substituted  cubes 
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have  lower  strain  energies.  On  this  scale,  the  PM3  values  show  large  variations  with 
respect  to  all  other  methods. 

Figure  4-5  presents  the  AEisoctesTniC  trends  relative  to  CCSD(T).  Again,  the 
BLYP  and  B3LYP  results  show  large  variations  over  the  series,  while  the  MBPT(2) 
and  CCSD  results  are  more  consistent.  The  range  of  the  BLYP  and  B3LYP  values 
are  27.4  and  33,  respectively,  while  the  MBPT(2)  and  CCSD  methods  have  a  much 
smaller  range  of  4.5  and  7  kcal  mol-1,  respectively.  The  largest  difference  for  the 
DFT  methods  is  for  N8,  where  B3LYP  differs  from  the  CCSD(T)  result  by  25.4  kcal 
mol-1. 

Table  4-6.  Reaction  coefficients  for  the  isodesmic  reaction:  nitroazacubane  +  a 
CH4  +  p  NH3  ->  7  N2H4  +  8  CH3N02  +  e  C2H6  +  C  N(CH3)3 


nitroazacubane 

a 

0 

7 

S 

e 

c 

(CN02)8 

24 

0 

0 

8 

12 

0 

1-aza 

21 

0 

0 

7 

9 

1 

1,3-diaza  (trans) 

18 

4 

I 

3 

1 

G 

7 

4 

I 

1,8-diaza  (opp) 

18 

1 

6 

7 

1,2-diaza  (cis) 

18 

4 
3 

1 

6 

7 

3 

1,3,5-triaza 

15 

0 

0 

5 

3 

3 

1,2,5-triaza 

15 

4 

1 

5 

4 

7 
3 

1,2,3-triaza 

15 

3 

2 

5 

5 

5 
3 

1,3,5,7-tetraaza 

12 

0 

0 

4 

0 

4 

1,2,3,5-tetraaza 

12 

8 
3 

2 

4 

2 

8 
3 

1,2,5,8-tetraaza 

12 

4 

3 

4 

3 

2 

1,2,3,7-tetraaza 

12 

8 
3 

2 

4 

2 

8 
3 

1,2,3,6-tetraaza 

12 

4 

3 

4 

3 

2 

1,2,3,4-tetraaza 

12 

16 
3 

4 

4 

4 

4 
3 

1,2,3,5,7-pentaaza 

9 

1 

3 

3 

0 

3 

1,2,3,5,6-pentaaza 

!) 

1C 
3 

4 

3 

1 

7 

1,2,3,4,5-pentaaza 

9 

20 
3 

5 

3 

2 

5 

1,2,3,4,5,7-hexaaza 

(trans) 

6 

8 

G 

2 

0 

2 

1,2,3,5,6,8-hexaaza 

(opp) 

6 

8 

G 

2 

0 

2 

1,2,3,4,5,6-hexaaza 

(cis) 

6 

28 

3 

7 

2 

1 

4 

Q 

1,2,3,4,5,6,7-septaaza 

3 

12 

9 

1 

0 

O 
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N8 

0 

16 

12 
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For  the  nitroazacubanes,  the  heats  of  formation  are  presented  in  Table  4-7  and 
graphically  in  Figure  4  6.  Beginning  with  octanitrocubane,  AH/s  progressively 
increase  along  the  series  with  increasing  nitrogen  content  of  the  cube.  When 
compared  to  the  azacubanes,  the  difference  in  the  MBPT(2)  and  SCF  values  is 
much  larger  due  to  the  presence  of  more  nitrogen  atoms  with  lone  electron  pairs. 
Differences  in  AH/  relative  to  MBPT(2)  are  presented  in  Figure  4-7.  All  three  of 
the  models  (SCF,  BLYP,  and  B3LYP)  display  the  same  trend  with  large  variations 
over  the  series.  One  interesting  aspect  in  the  energies  is  the  orientation  and 
repulsion  of  NO2  groups.  All  structures  have  been  fully  optimized  and  have  no 
imaginary  frequencies  at  the  SCF  level  of  theory.  Although  there  is  a  low  barrier  for 
rotation  of  a  N02  group,  the  repulsion  of  adjacent  nitro  groups  would  involve  a 
much  higher  energy.  For  1,2,3,4-tetraaza,  there  is  a  large  uncharacteristic  difference 
in  the  methods  with  respect  to  MBPT(2),  apparently  due  to  the  difficulty  in 
describing  the  repulsion  of  adjacent  nitro  groups. 

The  isodesmic  reaction  energies  for  the  nitroazacubanes  are  presented  in  Table 
4-8  and  graphically  in  Figure  4-8.  Again,  -Eisodesmic  decreases  along  the  series  with 
increasing  nitrogen  content  of  the  cube.  In  fact,  the  Eisodesmic  values  for  the 
nitroazacubanes  are  very  close  to  those  for  the  azacubanes.  Another  aspect  of  the 
nitroazacubane  energies,  however  is  the  orientation  and  repulsion  of  the  -N02 
groups.  There  is  a  small  peak  in  the  isodesmic  reaction  energies  for  the  molecule 
1,2,3,7-tetraaza.  Although  not  all  of  the  nitro  groups  are  attached  to  the  same  face 
of  the  cube  (see  Figure  4-1),  this  structure  seems  to  have  a  large  amount  of 
repulsion. 

Heats  of  formation  were  also  computed  directly  from  the  standard  states  of  the 
elements.  The  values  for  the  azacubanes  are  presented  in  Table  4-9  and  those  for 
the  nitroazacubanes  are  presented  in  Table  4-10.  The  differences  of  the  direct 
values  with  respect  to  the  most  accurate  isodesmic  method  are  presented  in  Figures 
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Table  4-9.  Heats  of  formation  in  kcal  mol  1  for  the  nitroazacubanes  computed 
directly. 


azacubane 

PM3 

BLYP 

B3LYP 

MBPT(2) 

(CN02)8 

205.1 

162.7 

181.0 

181.7 

1-aza 

203.4 

171.6 

187.4 

193.4 

1,3-diaza  (trans) 

204.0 

179.9 

193.0 

203.7 

1,8-diaza  (opp) 

202.0 

182.0 

195.6 

206.8 

1,2-diaza  (cis) 

215.1 

194.4 

208.5 

220.9 

1,3,5-triaza 

206.8 

187.7 

197.8 

212.7 

1,2,5-triaza 

216.0 

204.0 

215.5 

232.7 

1,2,3-triaza 

227.6 

215.9 

228.0 

246.3 

1,3,5,7-tetraaza 

212.2 

194.9 

201.8 

220.9 

1,2,3,5-tetraaza 

231.1 

224.9 

234.2 

256.8 

1,2,5,8-tetraaza 

230.7 

227.3 

237.0 

260.1 

1,2,3,7-tetraaza 

241.1 

236.4 

246.4 

270.0 

1,2,3,6-tetraaza 

240.7 

238.8 

249.2 

273.4 

1,2,3,4-tetraaza 

251.0 

251.6 

262.5 

288.1 

1,2,3,5,7-pentaaza 

247.3 

244.9 

252.0 

279.4 

1,2,3,5,6-pentaaza 

256.5 

261.1 

269.4 

299.1 

1 ,2,3,4,5-pentaaza 

265.0 

273.4 

282.3 

313.4 

1,2,3,4,5,7-hexaaza  (trans) 

282.3 

294.6 

301.4 

337.3 

1,2,3,5,6,8-hexaaza  (opp) 

283.4 

296.4 

303.5 

339.7 

1,2,3,4,5,6-hexaaza  (cis) 

290.3 

309.4 

317.1 

354.9 

1 ,2,3,4,5,6,7-septaaza 

317.3 

344.3 

350.7 

394.6 

N8 

352.8 

393.6 

399.3 

451.0 

4-9  and  4-10  for  the  azacubanes  and  nitroazacubanes,  respectively.  All  methods 
vary  over  a  large  range  and  the  accurate  description  of  AH;  directly  seems  to  be 
beyond  the  current  limits  of  ab  initio  correlated  methods.  Thus,  calculation  of  heats 
of  formation  directly  from  the  elements  is  not  the  preferable  route. 

4.4     Conclusions 
Although  a  caged  structure  might  be  difficult  to  achieve  synthetically, 
molecules  of  this  type  have  some  of  the  highest  densities.  Octanitrocubane's  density 
has  been  estimated  as  1.985  g  cm-3  [143].  All  of  the  azacubanes  and  the 
nitroazacubanes  are  highly  energetic  molecules.  With  increasing  nitrogen  content, 
greater  heats  of  formation  are  achieved  as  well  as  more  stable  structures  with 
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Table  4-10.  -AEisodesmic  in  kcal  mol  *  for  the  nitroazacubanes  computed  directly. 


nitroazacubane 

PM3 

BLYP 

B3LYP 

MBPT(2) 

(CN02)8 

172.2 

153.5 

207.9 

168.0 

1-aza 

178.4 

163.0 

210.2 

184.0 

l,3-diaza(trans) 

188.2 

173.7 

213.7 

199.5 

l,8-diaza(opp) 

183.0 

173.7 

214.1 

201.2 

l,2-diaza(cis) 

196.1 

189.4 

230.8 

218.4 

1,3,5-triaza 

200.0 

184.5 

216.9 

213.1 

1,2,5-triaza 

204.3 

199.7 

234.0 

233.0 

1,2,3-triaza 

215.6 

216.8 

248.8 

248.2 

1,3,5,7-tetraaza 

215.6 

193.7 

218.3 

223.9 

1,2,3,5-tetraaza 

229.3 

223.6 

251.4 

260.9 

1,2,5,8-tetraaza 

226.4 

226.6 

254.9 

265.2 

1,2,3,7-tetraaza 

233.3 

237.3 

266.5 

277.9 

1,2,3,6-tetraaza 

237.2 

235.9 

264.5 

274.8 

1,2,3,4-tetraaza 

244.1 

251.6 

281.5 

293.7 

1,2,3,5,7-pentaaza 

252.3 

244.9 

265.9 

284.3 

1,2,3,5,6-pentaaza 

258.0 

262.3 

285.0 

306.1 

1,2,3,4,5-pentaaza 

262.6 

273.8 

297.2 

319.4 

l,2,3,4,5,7-hexaaza(trans) 

287.8 

295.6 

312.3 

342.9 

l,2,3,5,6,8-hexaaza(opp) 

288.9 

299.6 

316.7 

348.3 

l,2,3,4,5,6-hexaaza(cis) 

291.7 

310.7 

328.4 

361.3 

1,2,3,4,5,6,7-septaaza 

321.3 

345.5 

357.0 

398.8 

Ns 

352.8 

393.7 

399.3 

450.9 

respect  to  fragments  of  the  isodesmic  reaction.  Semi-empirical  methods  are 
inadequate  for  this  problem,  though  PM3  is  better  than  the  others. 
Reparameterized  semi-empirical  theory  as  manifested  in  the  transfer  Hamiltonian 
might  be  able  to  overcome  these  limitations  [144]. 

Although  N8  is  too  unstable  for  production  in  bulk  quantities,  molecules  within 
the  azacubane  and  nitroazacubane  series  might  offer  more  kinetic  stability,  albeit  a 
lower  heat  of  formation.  Some  estimates  of  decomposition  pathways  and  energy 
barriers  for  the  azacubanes  have  been  made  [15,145]. 
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CHAPTER  5 
CHOICES  FOR  THE  ORBITAL  SPACE 

One  of  the  major  limitations  for  quantum  chemical  calculations  is  the 
dimension  of  the  molecular  orbital  basis.  In  fact,  the  computational  cost  for  ab 
initio  correlated  methods  has  a  tremendous  dependence  on  the  size  of  the  molecular 
orbital  basis.  Despite  many  chemically-interesting  processes,  it  is  prohibitively 
expensive  to  describe  systems  with  more  than  about  10  first- row  atoms  or  about 
300  molecular  orbitals  using  highly  correlated  methods. 

Typically,  molecular  orbitals  are  defined  from  canonical  Hartree-Fock  (HF) 
theory.  Although,  the  HF  determinant  has  the  lowest  energy  as  an  expectation  value 
of  a  single  Slater  determinant  with  the  Hamiltonian,  the  excited  or  'virtual'  orbitals 
are  purely  a  by-product  of  the  HF  calculation  in  a  basis  set,a   so  they  are  not 
necessarily  the  best  set  for  correlated  calculations.  In  this  chapter,  we  review  the 
merits  of  Hartree-Fock  theory,  and  provide  alternate  choices  for  molecular  orbitals. 

5.1     SCF  Orbitals 

In  Hartree-Fock  theory  [146],  the  exact  Hamiltonian  is  approximated  by  a  sum 
of  one-electron  Fock  operators 

7  7  n 

H  =  E  MO  +  E  r^  +  E  if1  ~  E  ^(»)  +  Nuclear  repulsion  (5-1) 

i  i<j  a</3        a@  i 

T(i)  =  h{z)  +  veff(i)  (5-2) 


a  For  atomic  and  diatomic  systems,  it  is  possible  to  solve  the  Hartree-Fock  equa- 
tions numerically,  rather  than  in  a  finite  set  of  Slater-  or  Gaussian-type  functions.  In 
these  numerical  solutions,  virtual  orbitals  are  not  defined. 
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The  set  of  orbitals,  {<pn},  are  optimized  to  give  the  lowest  electronic  energy  while 
remaining  orthonormal.  Such  conditions  are  imposed  using  Lagrange  multipliers 
and  functional  derivatives. 

c  [{<pn}]  =  e0  [{<?„}]  -  y,  A^  «»'W  -  <y  (5-3) 

Since  C  is  real  and  {i\j}*  =  (j\i),  the  multipliers  constitute  an  Hermitian  matrix. 

AJ  =  Ay  (5-4) 

Taking  the  partial  derivative  of  C  with  respect  to  ip*k  yields: 
—  =  (8<pk\h\<pk)  +  -  ^2{6(pk<pj\\<Pk'Pj)  +  2  Yl{<Pis<Pk\\<Pi<Pk)  -^2hJ(S<Pk\<Pj) 

J  i  j 

(5-5) 
Setting  Equation  5-5  equal  to  zero,  factoring  8ipk  from  every  term,  and  combining 
the  summation  indices  gives: 


h\(pk)  + 


E  (*  -  *>) 


l^)-EA^l^)=°  (5-6) 


where  J  and  K  have  been  introduced  to  denote  the  Coulomb  and  exchange 
operators,  respectively.  Because  ve^  =  V  •  ( Jj  —  Kj  1 ,  Equation  5-6  is  a 
pseudo-eigenvalue  equation  in  terms  of  the  Fock  operator.  It  is  often  beneficial  to 
diagonalize  the  A  matrix  to  provide  an  energy  for  each  orbital. 

e  =  U^XU  (5-7) 

However,  there  are  several  considerations  to  obtaining  orbital  energies.  The 
first  is  whether  expectation  values  will  differ  between  the  two  wave  functions.  If  A  is 
diagonalized  by  a  real  unitary  matrix  (Uf  =  U_1;  L^U  =  l),  the  new  wave  function 
can  be  written  in  terms  of  the  old  one  as 

|tf')  =  de*(U)|tt)  (5-8) 


G5 
and  the  expectation  values  as 

<tf'|i|^'>  =  det  (UfU)  (#|i|#)  =  det  (1)  (tf|i|tf )  (5-9) 

Hence,  all  expectation  values  will  be  invariant  to  unitary  transformations  of  the 
wave  function.  The  second  question  is  whether  the  Fock  operator  is  invariant.  The 
only  part  of  the  Fock  operator  that  depends  on  the  orbitals  is  ve^  and  it  is 
considered  below: 

^'(l)     =    E/^(2)~%i(2)^  (5-10) 

3 

=  E/E^^^E^2^ 

j  J     k  12       I 

=    51  [  Swpffl— —<Pi(2)dr2 

ki  J  ri2 

=    Y,  l^2)1-— —M2)dr2 
kJ  ri2 

Consequently,  we  can  write  the  SCF  canonical  equation  as 

:F(1)¥>„(1)  =  ev<pv(l)  (5-11) 

However,  thirdly,  the  canonical  orbitals  are  not  typically  localized  in  a  region  of 
space,  but  rather  spread  over  many  atoms.  This  complicates  the  interpretation  of 
molecular  orbitals  in  terms  of  chemical  bonds  which  are  mostly  localized. 

One  of  the  features  of  the  canonical  orbital  energies  is  that  they  may  be 
interpreted  as  approximate  ionization  potentials.  This  is  Koopmans'  theorem  and  it 
involves  the  same  unrelaxed  molecular  orbitals  for  the  n-1  system  as  in  the  n 
electron  system.  For  instance,  the  energy  of  an  n  electron  system  is 

n  n 

E{n)  =  52(^1*)  +  2  5>illV>  (5-12) 


GG 
Provided  the  orbitals  are  frozen,  the  energy  is  very  similar  for  the  n-1  system. 

E(n-l)  =  ^2(i\h\i)  +  lf^(ij\\ij)  (5-13) 

i  '    i,j 

The  ionization  potential  (IP)  is  the  difference  between  the  n-1  and  n  electron 
energies 

IP  =  E(n)  -  E(n  -l)  =  -hn-U  Y^(nj\\nj)  +  ^HH  )  =  ~e«       (5-14) 

which  is  also  the  orbital  energy  for  whichever  electron  is  removed. 

5.2     Vn-i  and  Vn-a  Virtual  Orbitals 
With  some  inspection,  it  is  easy  to  realize  that  the  character  of  the  occupied 
and  virtual  elements  of  the  Fock  matrix  are  quite  different. 

n 

Fpq  =  hpq  +  ^2(ip\  \iq)  (5-15) 

i 

For  the  occupied  orbitals,  the  effective  potential  terms  cancel,  so  they  are 
determined  in  the  potential  of  n-1  electrons.  However,  for  the  virtual  orbitals,  p  and 
q  are  never  equal  to  i,  so  they  are  determined  in  the  potential  of  n  electrons,  and 
hence,  might  be  more  appropriate  for  the  n+1  electron  system.  Furthermore,  since 
the  HF  energy  is  invariant  to  any  definition  of  the  virtuals,  it  might  be  beneficial  to 
redefine  them  [147].  For  CC  or  CI  methods  that  include  all  excitations  of  single, 
double,  etc.  type,  the  computed  energies  and  properties  will  be  invariant  to  virtual 
orbital  transformations.  However,  the  invariance  is  lost  if  only  a  subset  of  the 
virtuals  are  retained.  Thus,  it  is  possible  to  define  better  virtuals  where  subsets 
retain  most  of  the  correlation  energy,  as  they  require  much  less  time  than  the 
calculations  in  the  full  virtual  space. 

Kelly  [148-152]  investigated  excluding  different  orbitals  from  ve//  to  determine 
virtual  orbitals  and  MBPT  correlation  energies  for  various  orders.  It  is,  however, 
somewhat  ambiguous  which  occupied  orbital  should  be  excluded.  A  more 
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satisfactory  route  in  modifying  ve^  is  to  exclude  each  occupied  orbital  in  an 
averaged  way  [153]. 

K-i  =  —v*ff  (5-16) 

It  is  also  possible  to  envision  excluding  a  larger  number  of  electrons  or  even 
fractional  numbers  of  electrons.  There  ultimate  justification  will  be  that  a  larger 
amount  of  correlation  energy  is  obtained  in  a  smaller  subspace.  If  this  is  done  in  the 
same  averaged  way,  the  new  potential  has  the  form 

Vn-a  = veff  (5-17) 

Numerical  results  for  the  Vn_x  and  Vn_a  potentials  where  a  is  the  number  of  valence 
electrons  are  presented  in  Section  5.5. 

5.3     Density  Functional  Theory  Virtual  Orbitals 
Although  originally  formulated  for  solids,  today,  density  functional  theory 
(DFT)  is  routinely  applied  to  atoms  and  molecules  [154].  The  one-particle 
Kohn-Sham  equations  are  very  similar  in  form  to  the  Hartree-Fock  equations, 
except  that  the  Kohn-Sham  effective  potential  is 

veff(r)  _  9J[p]      dExe\p] 

v  {r)  -  dP(r) +  ~d^y  (   8) 

In  Kohn-Sham  theory,  ve//  contains  Coulomb  repulsion,  exchange,  and  correlation, 
and  in  principle  is  exact  provided  the  exchange-correlation  functional  is  correct. 
One  open  question  is  how  DFT  virtual  orbitals  will  perform  in  high-level  CC 
calculations,  since  unlike  SCF  orbitals  an  electron  in  an  occupied  or  virtual  orbital 
feels  the  same  potential.  In  Section  5.5,  CC  correlation  energies  using  virtual 
orbitals  from  the  LDA/VWN  [155, 156]  and  P W91  [157]  functionals  are  presented 
and  compared  with  other  choices  for  the  virtual  space. 
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5.4     Frozen  Natural  Qrbitals 

Natural  orbitals  were  introduced  by  Dirac  [158]  and  extensively  studied  by 
Lowdin  [159].  They  are  denned  as  the  eigenvectors  of  the  one-particle  density 
matrix  and  provide  the  fastest  convergence  of  the  configuration-interaction  (CI) 
expansion  [159,160].  One  of  the  potential  problems  with  natural  orbitals  is  that  a 
correlated  wave  function  is  needed  for  their  determination.  (They  can  also  be 
determined  from  orthogonalizing  the  Feynman-Dyson  amplitudes  using  some 
self-energy  approximation  [161,162]).  Several  authors,  including  Meyer  [163],  used 
pseudo-natural  orbitals  obtained  from  an  approximate  density  matrix  at  a  lower 
level  of  theory  to  generate  the  natural  orbitals.  Since  the  density  matrix  is  from  a 
lower  level  of  theory,  the  generation  of  pseudo-natural  orbitals  does  not  involve  any 
rate  limiting  steps  or  is  not  a  bottleneck  for  the  calculation. 

In  practice,  it  is  advantageous  if  there  are  no  off-diagonal  occupied-occupied, 
fjj,  or  virtual-virtual,  iab,  elements  in  the  Fock  matrix,  since  these  can  force  an 
iterative  solution.  For  example  with  CCSD(T),  the  T3  amplitudes  are  given  by  [164] 

D$£$k    =    Y,P(i/jk)P(a/bc)t%(bc\\ei)  -  J3P(i/i*)P(o/6c)l*S,0'*l|mo> 

e  m 

+P(i/jk)P(a/bc)  [tUbc\\jk)  +  fiat%] 

-P(i/jk)  x  £(l  -  8im)ftmt%  x  £(1  -  8ae)faet\%  (5-19) 

m  e 

The  last  terms  in  Equation  5-19  contain  fim  and  fae  elements  coupled  to  T3.  If  the 
fim  or  fae  elements  are  not  zero,  the  T3  amplitudes  appears  on  both  sides  of  the 
equation  and  it  must  be  solved  iteratively.  However,  we  can  use  the  invariance  of 
CC  theory  to  orbital  rotations  among  the  occupied  and  virtual  spaces  to  choose  to 
make  a  semicanonical  transformation  of  the  resultant  Fock  matrix  to  simplify  the 
solution  of  Equation  5-19.  With  semicanonical  orbitals,  the  ftj  (i^j)  and  fab  (a^b) 
elements  are  zero,  however,  the  Fock  matrix  does  contain  fia  terms.  Since  the  fia 
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elements  are  coupled  to  T2  amplitudes  instead  of  T3's  in  Equation  5-19,  they  do 
not  force  an  iterative  solution. 

One  close  approximation  to  natural  orbitals  that  avoids  the  complications  of  % 
(i#j).  Ub  (a^b),  and  even  fio  terms  in  energy  calculations  is  frozen  natural  orbitals 
(FNOs)  [165].  In  this  approach,  only  the  virtual-virtual  block  of  the  one  particle 
density  matrix  is  diagonalized  to  generate  natural  orbitals  constrained  to  the  virtual 
space.  Within  the  full  virtual  space,  FNOs  are  just  a  unitary  transformation  of  the 
original  (e.g.,  Hartree-Fock)  orbitals,  so  the  energy  is  invariant.  However,  FNOs  can 
achieve  nearly  all  of  the  energy  in  a  smaller  space  which  provides  considerable 
savings  in  the  resulting  calculation.  Since  some  virtual  orbitals  are  excluded  from 
the  resulting  calculation  or  dropped,  the  energy  is  no  longer  invariant.  Furthermore, 
the  active  block  of  the  virtual-virtual  Fock  matrix  can  be  diagonalized,  so  there  are 
no  fo6  (a^b)  terms  which  need  to  included  in  the  energy  calculation.  However,  the 
Fock  matrix  is  not  block  diagonal  with  respect  to  the  active  and  inactive  portions, 
so  there  are  fa6  (a^b)  terms  within  the  inactive  space  which  are  pertinent  in 
gradient  calculations.  Some  numerical  results  of  FNOs  taken  from  a  MBPT(2)  wave 
function  are  shown  in  the  next  section. 

5.5     Illustrative  Examples 

Table  5-1.  Percent  of  total  CCSD  and  CCSD(T)  correlation  energies  recovered  for 
different  orbital  choices.  The  molecular  system  is  CO  with  r  = 
1.12832  A  [5]. 


46.6 


Percent  of  virtual  space  (number  of  virtual  orbitals) 

CCSD  CCSD(T) 

80%(82)          60%(62)  40%(4iy        80%           60%         ~40%" 

SCF                         89.4                77.0  46.7              89.2          76.4 

Vn-i                        90.0                78.9  53.6 

Vn-a                        95.7                90.1  48.8 

LDA/VWN             89.4                77.1  46.8 

PW91                      89.4                77.1  46.7 

FNQ       1000               99-4  97.1             100.0         99.4           96.9 


70 


100 


Bfl 

s- 

c 

d 

o 

s 
b 

o 
o 

(+- 
o 


90 


80 


70 


°    60 


d 

B 

u 

a. 


50 


40 


, 

x^Z^^^ 

-•-FNO 
-*-SCF 

LDA/VWN 
*  PW91 
-*-Vn-l 
-•-Vn-a 

yy 

y 

40 


50  60  70 

percent  of  virtual  space 


80 


Figure  5-1.  Percent  of  total  CCSD  correlation  energies  recovered  for  different 
orbital  choices.  The  molecular  system  is  CO. 

In  this  section,  a  number  of  examples  are  considered  to  show  the  numerical 
performance  of  different  orbital  choices  in  a  reduced  space.  In  Table  5-1,  the 
percentage  of  the  correlation  energy  is  tabulated  for  different  choices  of  the  orbital 
space  and  truncations  of  that  space.  The  molecule  is  CO  described  in  an 
aug-cc-pVTZ  basis  with  all  core  electrons  frozen.  The  results  are  also  presented 
graphically  in  Figure  5-1.  The  SCF,  LDA/VWN,  and  PW91  virtual  orbitals  recover 
nearly  identical  amounts  of  the  CCSD  correlation  energy  for  the  three  truncations 
of  the  space.  The  Vrn_1  and  Vn-a  results  recover  slightly  more  of  correlation  energy 
in  the  reduced  spaces.  However,  the  FNOs  perform  much  better  for  both  the  CCSD 
and  CCSD(T)  correlation  methods.  For  CCSD,  FNOs  achieve  100.0%,  99.4%,  and 
97.1%  of  the  correlation  energy  in  80%,  60%,  and  40%  of  the  virtual  space, 
respectively. 
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Table  5-2.  Percent  of  total  CCSD  and  CCSD(T)  correlation  energies  recovered  for 
different  orbital  choices.  The  molecular  system  is  C2H4  with  rcc  = 
1.3342  A,  vch  =  10812  A,  and  &Cch  =  121.3  [5]. 

Percent  of  virtual  space  (number  of  virtual  orbitals) 


CCSD 


CCSD(T) 


80%(162)         60%(121)  40%(81)  80%  60%         40% 


SCF 

95.6 

K-i 

95.8 

'  n— a 

97.3 

LDA/VWN 

95.5 

PW91 

95.5 

FNO 

100.0 

78.5 
79.9 
90.1 
77.0 
76.9 
99.8 


55.6 
61.6 
48.6 
56.0 
56.0 
98.6 


95.5 


77.8 


100.0 


99.8 


54.6 


98.5 


Table  5-3.  Percent  of  total  CCSD  and  CCSD(T)  correlation  energies  recovered  for 
different  orbital  choices.  The  molecular  system  is  F2  with  r  = 
1.41193  A  [5]. 

Percent  of  virtual  space  (number  of  virtual  orbitals) 

CCSD(T) 


CCSD 


80%(81)         60%(61)  40%(40)         80% 


60% 


40% 


SCF 

89.1 

K-i 

89.6 

Vn—a 

96.8 

LDA/VWN 

89.1 

PW91 

89.0 

FNO 

99.8 

74.6 
76.1 
90.0 
74.6 
74.5 
98.7 


52.8 
52.8 
53.2 
52.8 
52.8 
93.2 


88.9 


74.2 


99.8 


98.6 


52.5 


92.8 


Tables  5-2  through  5-5  present  additional  results  for  the  C2H4,  F2,  H20,  and 
NH3  molecular  systems.  One  additional  trend  is  that  the  Vn-a  potential  performs 
better  for  80%  and  60%  reductions  of  the  virtual  space,  but  the  Vn_i  potential 
performs  slightly  better  for  the  40%  reduction  of  the  virtual  space. 

Because  total  energies  are  seldom  important,  it  is  more  constructive  to  examine 
computed  energy  differences.  In  Table  5-6,  activation  energies  for  the  unimolecular 
dissociation  of  the  pentazole  anion,  N5,  are  presented.  Again  the  SCF,  LDA/VWN, 
and  PW91  results  are  quite  similar.  The  results  for  the  FNOs  are  remarkably  good. 
For  the  CCSD(T)  method,  FNOs  with  only  40%  of  the  virtual  space  predict  an 
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Table  5-4.  Percent  of  total  CCSD  and  CCSD(T)  correlation  energies  recovered  for 
different  orbital  choices.  The  molecular  system  is  H20  with  voh  — 
0.9572  A  and  aHOH  =  104.52  [5]. 

Percent  of  virtual  space  (number  of  virtual  orbitals) 


CCSD 


CCSD(T) 


80%  (80) 

60%  (60) 

40%  (40) 

SCF 

85.5 

73.8 

38.0 

Vn-l 

87.0 

77.1 

48.0 

'n— a 

95.8 

90.6 

4.7 

LDA/VWN 

85.5 

73.8 

38.0 

PW91 

85.5 

73.7 

38.2 

FNO 

100.0 

99.7 

98.3 

80% 


60% 


40% 


85.2 


73.2 


38.1 


100.0 


99.7 


98.2 


Table  5-5.  Percent  of  total  CCSD  and  CCSD(T)  correlation  energies  recovered  for 
different  orbital  choices.  The  molecular  system  is  NH3  with  rNH  = 
1.0116  A  and  &NHN  =  106.7  [5]. 

Percent  of  virtual  space  (number  of  virtual  orbitals) 


CCSD 


80%(100) 

60%(75) 

40%(50) 

SCF 

90.9 

77.9 

47.1 

K»-i 

93.0 

80.5 

59.0 

*n— a 

97.4 

92.2 

17.1 

LDA/VWI 

*             90.9 

77.8 

47.2 

PW91 

90.9 

77.8 

47.0 

FNO 

100.0 

99.8 

98.7 

CCSD(T) 


80% 


60% 


40% 


90.6 


100.0 


77.1 


99.8 


46.4 


98.6 


activation  energy  of  29.3  kcal  mol"1.  This  is  within  1  kcal  mol-1  of  the  CCSD(T) 
result  for  the  full  virtual  space,  28.8  kcal  mol-1. 
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CHAPTER  6 
CONCLUSIONS 

Theoretical  chemistry  methods  have  a  great  potential  to  study,  predict,  and 
discover  high-energy  density  materials.  Often  an  initial  theoretical  survey  can 
identify  several  promising  highly  energetic  molecules  that  may  be  investigated 
further.  In  the  more  detailed  work,  theory  can  closely  assist  the  experimental 
efforts.  In  fact  CC  methods  have  produced  accurate  structures  and  energetics, 
however,  one  of  the  bottlenecks  in  their  conventional  formulation  and  application  is 
the  dimension  of  the  molecular  orbital  basis. 

This  dissertation  has  explored  several  different  choices  for  the  orbital  space  and 
the  numerical  results  show  that  CC  methods  are  largely  invariant  to  different 
choices  for  the  virtual  orbital  space  including  DFT,  Vn-i,  and  Vn-a  potentials. 
However,  approximate  frozen  natural  orbitals  define  a  better  virtual  space  of 
reduced  dimension  where  much  more  of  the  correlation  energy  is  recovered.  The 
FNO  procedure,  detailed  in  Appendix  A,  constructs  a  reduced  space  with  natural 
orbital  information  from  the  full  space.  Hence,  FNOs  achieve  a  large  savings  in  the 
time  needed  per  iteration  of  the  CC  equations.  Although  not  investigated  in  this 
dissertation,  FNOs  might  be  useful  in  coupled-cluster  equation-of-motion 
calculations  for  electronically  excited,  ionized  and  attached  states.  Additional 
future  work  could  extend  analytical  gradient  techniques  for  references  with  FNO 
reductions  of  the  virtual  space. 
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APPENDIX  A 
COMPUTATIONAL  IMPLEMENTATION 

The  ability  to  generate  frozen  natural  orbitals  has  been  implemented  into  the 

ACES  II  suite  of  computer  programs  [68].  The  following  is  a  summary  of  the 

procedures  performed  by  the  xfno  ACES  member  executable.  It  is  important  to 

note  that  in  the  limit  of  the  full  virtual  space  (i.e.,  5=0),  this  procedure  results  in 

canonical  Hartree-Fock  virtual  orbitals. 

•  The  virtual-virtual  block  of  the  relaxed  one-particle  density  matrix  is 
calculated.  For  MBPT(2),  this  is: 


M»6  —  7;  7  ,  7 — ; w — ; 7  (A_1) 


•  Frozen  natural  orbitals,  U,  and  their  occupation  numbers  are  obtained  by 
diagonalizing  this  matrix: 

U]DabU  =  u  (A-2) 

•  The  natural  orbitals  are  transformed  to  the  basis  of  symmetry-adapted 
orbitals. 

vrt        vrt         vrt 
so  C  vrt  U  =  so  U'  (A-3) 

•  The  virtual-virtual  block  of  the  Fock  matrix  is  built  in  the  basis  of  natural 

orbitals. 

s?+       so      vrt  vrt 

vrt  U  f  so  F  so  U  =  vrt  F'  (A-4) 

•  The  active  block  of  the  Fock  matrix  is  diagonalized  to  produce  new  orbital 
energies  for  the  reduced  space.  The  orbital  energies  for  the  reduced  space  are 
always  larger  than  the  original  orbital  energies. 

vrt  -6  vrt  -  S  vrt  _  g 

vrt-  8      Z*      vrt  -8      F'      vrt  -  8      Z     =  e'  (A-5) 

•  The  orbitals  are  updated  to  diagonalize  the  active  virtual  block  of  the  Fock 
matrix. 

vrt  -  8  vrt  -  5         vrt  -  8 

so      U'     vrt -8      Z     =so     U"  (A-6) 
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